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1.  INTRODUCTION 

This  report  presents  technical  matter  related  to  the  research  performed 
at  the  Department  of  Physics  of  the  University  of  Modena  for  the  Contract 
Number  DAJA45-86-C-0004  "Monte  Carlo  Analysis  of  Quantum  Transport  and 
Fluctuations  in  Semiconductors”.  Such  a  Contract  follows  a  similar  Contract 
with  the  same  title  relative  to  the  period  1983/85. 

In  the  Final  Technical  Report  of  that  first  Contract  we  reported  a  general 
review  of  quantum  approaches  to  charge  transport  in  semiconductors;  we  there¬ 
fore  refer  to  that  review  for  general  references  and  for  introductory  treatment 
of  quantum  transport. 

In  the  present  report  we  shall  limit  ourselves  to  the  results  of  the  research 
performed  under  the  Contract  to  which  it  relates.  For  the  sake  of  completeness, 
however,  v/e  shall  sometimes  include  some  materials  developed  during  the  first 
part  of  this  research. 

In  the  presentation  we  shall  follow  the  same  order  as  in  the  Work  Statement 
of  the  Contract.  In  particular,  Chapters  2,  3,  and  4  refer  to  the  points  1,  2,  and 
3  listed  as  Objectives  of  the  Research. 

As  it  regards  quantum  transport,  two  major  lines  of  activity  have  been 
followed.  In  one  of  them,  described  in  Sec.2.1,  the  traditional  Monte  Carlo 
approach  has  been  extended  to  include  the  quantum  features  that  could  be 
represented  by  the  use  of  the  joint  spectral  density  in  the  transition  probabilities. 
The  theory  at  the  basis  of  this  extension  has  been  developed  to  a  certain  level 
of  detail  in  general,  and  in  some  special  cases  it  has  been  fully  developed  up 
to  the  realization  of  Monte  Carlo  codes  for  simple-model  semiconductors(Secs. 
2.1.4,  2.1.5). 

In  the  second  approach  a  method  is  developed  for  the  numerical  solution 
of  the  Liouville-von  Neumann  equation  that  describes  electron  transport  in 
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semiconductors  in  a  fully  quantum  mechanical  formulation.  The  problem  has 
been  solved  in  principle,  but  the  amount  of  computations  necessary  to  obtain 
the  results  for  real  cases  of  interest  is  still  very  large.  The  method  is  presented  in 
Sec.2.2.  Two  applications  have  been  considered  in  particular:  the  analysis  of  the 
effect  of  multiple  collisions  and  of  intra-collisional  field  effect  in  presence  of  high 
applied  electric  fields  (Sec. 2. 2. 4.1),  and  the  analysis  of  the  energy  relaxation  of. 
photoexcited  electrons  in  absence  of  external  fields  (Sec. 2.2. 4. 2). 

A  comparison  of  the  two  approaches  presented  in  this  report  should  stress 
that  the  first  one  requires  little  modifications  of  the  traditional  Monte  Carlo 
technique, but  some  more  investigations  are  required  to  understand  how  much  of 
the  quantum  features  of  electron  transport  it  will  be  able  to  describe;  the  second 
approach  is  more  rigorous,  but  further  analytical  developments  are  required  if 
we  want  it  to  be  useful  for  practical  calculations. 

As  it  regards  the  analysis  of  correlation  functions  of  hot  electrons  (Chapter 
3),  two  major  lines  of  research  have  again  been  followed. 

In  the  first  work  (Sec. 3.1),  explicit  formulae  have  been  obtained  for  the 
time  evolution  of  the  velocity  correlation  functions  which  generalize  the  linear- 
response  theory  to  conditions  far  from  equilibrium. 

The  second  line  of  research  (Sec.  3.2)  has  investigated  the  velocity  corre¬ 
lation  function  of  electrons  in  semiconductor  quantum  wells  under  high  applied 
electric  fields  both  in  transient  and  in  steady-state  conditions. 

Hot  phonons,  i.e.  phonon  populations  in  excess  with  respect  to  their 
thermal-equilibrium  values,  can  be  obtained  by  means  of  the  application  of 
strong  electric  fields  or  by  photoexcited  high-energy  electrons  with  laser  pulses. 
This  subject  has  recently  collected  large  interest  and  debates.  Chapter  4  con¬ 
tains  the  description  of  the  method  developed  by  our  group  for  a  simultaneous 
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simulation  of  electron  and  phonon  transport  in  bulk  semiconductors  and  quan¬ 
tum  wells.  Several  applications  are  discussed  in  details  for  bulk  GaAs  and  for 
GaAs-AlGaAs  quantum  wells.  Good  qualitative  agreement  with  Raman  exper¬ 
imental  data  is  obtained. 


2.QUANTUM  TRANSPORT 


2.1  Inclusion  of  the  Spectral  Density  in  the  conventional  Monte  Carlo 

Approach 

2.1.1  Introduction 

The  fast  development  in  the  field  of  submicron  devices  has  provided  a 
renewed  interest  in  the  theory  of  electron  transport  beyond  the  free-particle 
approach  based  on  the  ser  [classical  Boltzmann  equation.  Indeed,  quantum 
theory  has  indicated  that  a  proper  treatment  of  high-field  transport  in  semi¬ 
conductors  should  include  the  intra  collisional  field  effect  (ICFE)  as  well  as 
collisional  broadening  (CB)  [2.1]. 

ICFE  accounts  for  the  presence  of  the  electric  field  in  the  collision  operator 
of  the  .kinetic  equation.  In  other  words,  a  scattering  event  does  not  occur 
between  states  described  by  the  plane  waves  of  a  free  electron,  but  between 
those  of  an  electron  in  the  field.  Introduced  by  Levinson  and  Yasevichyute 
[2.2j,  ICFE  has  subsequently  been  investigated  by  many  researchers  [2.3  -  11], 

Since  the  total  Hamiltonian  H  is  not  the  same  as  that  of  the  perfect  crystal, 
an  electron  in  a  state  with  wave  vector  k  is  not  in  an  eigenstate  of  H  and  does 
not  correspond  to  a  well  defined  energy  hui,  but  to  a  broad  interval  of  energies 
associated  to  k.  The  spectral  function  A(k,w)  gives  the  probability  that  an 
electron  in  state  k  is  found  with  energy  hu.  It  is  possible  to  show  that  A(k,w) 
represents  also  the  probability  that  an  electron  with  energy  hwis  found  is  state 
k  [2.12,13]. 

The  effective  mass  model  is  recovered  when  A(h,  w)  =  27t6(h2k2 /2m  -  hui).  CB 
in  the  context  of  high-field  transport  in  semiconductors  was  described  in  detail 
by  Barker  [2.14]  and  its  importance  within  a  Monte  Carlo  simulation  was  firstly 
pointed  out  by  Capasso  et  al  [2.15], 
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Since  then,  several  ’esearch  groups  [2.16  -  24]  have  attempted  to  esti¬ 
mate  the  importance  of  these  effects  by  suitable  generalizations  of  standard 
Monte  Carlo  algorithms.  However,  the  technical  details  have  been  scarcely 
documented.  As  a  consequence,  it  has  been  practically  impossible  for  anyone 
outside  a  given  research  group  to  reproduce  quantitatively  the  results  obtained 
by  other  groups.  This,  in  turn,  has  led  to  contradictory  findings  and  difficulties 
in  objectively  assessing  the  significance  of  the  various  quantum  effects,  and  in 
judging  the  relative  merits  of  various  calculations. 

The  main  objecrive  of  this  section  is  thus  to  present  the  theoretical  frame- 
wotk  and  the  detail-  cf  two  Monte  Carlo  algorithms  which  we  have  used  to  esti¬ 
mate  ICFE  and  CB  effects  within  a  simple  model  semiconductor.  Our  scheme  is 
based  on  a  quantum  kinetic  equation  we  have  recently  derived  [2.25],  and  which 
accounts  for  quantum  effects  through  the  introduction  of  the  joint  spectral  den¬ 
sity  K(ti,Cf)  describing  the  connection  between  the  initial,  e,-,  and  final,  £/, 
kinetic  energy  of  the  carrier  during  a  single  scattering  process.  (Let  us  remark 
that  we  are  within  the  so  called  completed  collision  regime  [2.11.25],  which 
enables  us  to  construct  time-independent  scattering  rates).  The  joint  spectra! 
density  reduces  to  the  •veil  known  delta  function  behavior  when  ICFE  or  CB 
are  neglected.  Therefore,  the  present,  scheme  yields  the  correct,  semiclassical 
Boltzmann  limit.  One  of  the  objective  of  our  study  is  to  provide  a  first  stan¬ 
dard  on  the  subject  which  hopefully  will  be  of  general  use  and  open  to  further 
improvements. 

2.1.2  Model 

The  theory  is  formulated  in  terms  of  the  Generalized  Kadanoff-Baym  ap¬ 
proach  [2.26,  27].  Our  physical  model  is  based  on  the  following  assumptions: 

(i)  One  isotropic  band; 

(ii)  Time  independent  a  id  space  homogeneous  conditions; 
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(iii)  An  ansatz  which  links  the  single  particle  Wigner  distribution  function  /( k) 

to  the  full  correlation  function  G<(k,t)  [2.28]  (G<(k,t)  =  r'A(k,t)/(k-  |t|)); 

(iv)  Completed  collisions  limit. 

Within  this  model,  the  quantum  kinetic  Boltzmann  equation  for  /(k)  writes 
[2.25]: 

/00  =  /  dt  [w«M(k'  -  eEt/ft,k  -  eEt/h)  f(k'  -  eEt/h) 

-  WQM( k  -  eEt/h,  k'  -  eEt/h) /(k  -  eEt/A)]  (2.1) 

The  scattering  probability  per  unit  time  W^M ,  for  electron  phonon  interactions, 
is  given  by: 

^gM(kx,k2)  =  lV^)|S  (N,  +  \  +  \v)  Jf(k!,ka)  (2.2) 

»7  =  ±1 

where  /f(ki,k2)  is  the  joint  spectral  density  of  the  quasi  particle  which,  in  the 
self-consistent  Born  approximation,  is  given  by: 

/°o  _-p  .p 

dt'2Re{A{  k!  +  — *',t'M(ka  +  —  t',-t')exp(-xr,u/)} 

(2.3) 

The  description  of  the  physical  processes  is  based  on  a  quasi-particle  picture, 
where  k  and  k'  are  the  quasi  particle  kinetic  momenta  before  and  after  a  scat¬ 
tering  event,  q  =  |k~k'|  is  the  transferred  momentum,  E  is  the  external  applied 
electric  field,  |k,(<z)|2  the  square  of  the  matrix  element  for  electron  phonon  col¬ 
lision  ,  Nq  the  equilibrium  phonon  population,  hwq  the  phonon  energy  involved 
in  the  collision,  with  r)  =  ±1  referring  respectively  to  emission  and  absorption 
processes. 

The  joint  spectral  density  /f(ki,k2)  is  the  central  quantity  in  our  ap¬ 
proach,  since  it  enables  us  to  account  for  quantum  corrections  of  the  otherwise 
8emiclassical  free-particle  picture. 
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2.1.3  The  Spectral  Density 

In  the  framework  of  a  many-body  approach  which  uses  a  the  Green  function 
formalism  [2.29],  the  spectral  density  is  defined  as: 

A(k,t)=»[G'(k,t)-G‘(k,t)]  (2.4) 

where  Gr,a(k,<)  is  the  retarded  (advanced)  Green  function,  respectively. 

In  the  presence  of  a  steady  external  electric  field  E  and  of  collisions,  Gr(k,t) 
obeys  the  Dyson’s  equation  which,  when  expressed  in  terms  of  gauge  invariant 
variables,  has  the  form: 

Gr(k,t)  =Gi(k,t)  +  J  J  dtldt2GrE{k-e^(t-tl),t1)x 

x  Er(k  -  —  (*2  -  fi),f  -  (ti  +  t2))Gr(k  -  —  («2  -  0>*2)  •  (2.5) 

Er(k,  t)  is  the  retarded  self  energy.  The  field  dependent  retarded  Green  function 
Gg(k,t)  is  given  by  [2.26]: 

rt/2 

Gg(k,<)  =  -i8(t)exp  [-f  /  du  t(k  -  eEu/fr)]  (2.6) 

J-t/2 

where  6(t]  is  the  unit  step  function,  and  e  is  the  carrier  kinetic  energy. 

An  equivalent  formulation  of  Eqs. (2.4-6)  in  the  frequency  domain  can  be 
easily  obtained  through  a  Fourier  transformation. 

2.1.4  Applications 

In  the  following,  the  explicit  expressions  for  Ff(k,k')  =  Ff(c,,  £/)  will  be 
given  for  the  case  of  ICFE  and  CB. 


2. 1.4.1  Free  Electrons  with  ICFF- 

For  free  electrons  in  the  presence  of  an  external  electric  field  E  the  joint 
spectral  density  is  given  by  (2.30]: 


♦*(£)  [■-"(«&)]} 


with 


P  =  £/  —  c,-  +  hu„ ;  Q  — - -eq  •  E, 

m 


where  C  =  C(z)  and  S  =  S[z)  are  the  Fresnel  integrals  given  by: 
C(z)  =  (-)5  f  dtcos(t2);  S(z)  =  (— )s  f  dt  sin(t2) . 

x  J  o  *  Jo 


(2.7) 


(2.8) 


(2.9) 


Because.of  the  presence  of  fast  oscillations  associated  to  the  Fresnel  integrals  in 
Eq.  (2.9)  A' («,-,£/)  in  the  form  of  Eq.(2.7)  cannot  be  used  directly  in  a  Monte 
Cj  rlo  scheme  which  requires  a  positive  definite  quantity.  As  suggested  by  Barker 
[2.3],  a  plausible  way  of  suppressing  the  oscillations  (whose  tails  integrate  to 
zero)  is  to  approximate  the  expression  (2.7)  with  a  Lorentzian  function.  In 
order  to  preserve  normalization  and  control  the  low  and  high  energy  tails  we 
use  a  truncation  procedure:  K(u, «/)  =  0  for  tj  <0  and  tj  >f/,  where  Ij  is  a 
cut-off  energy.  Thus,  Eq.  (2.7)  is  approximated  by  the  following  expression: 


*(«,«/) 


1 _ A 

|  0  \l/ 2  1  +  ( Z,  -XJ-XO-  7T|^)3 


(2.10) 


where 


A  =  [arctg(x/  -  n  +  x0  +  Jr—j) 


Q 

•  arctg{x0  -  x<  +  *'j-^-|)]~1 


(2.11) 


0  = 


chE 


jr(2m*)VJ 


( xlj2cos0f  -  xl'^cosOi) 


(2.12) 
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Here  6{  and  6 j  are  respectively  the  initial  and  final  angles  of  k,  and  k/  with 
respect  to  the  electric  field  E  and  we  introduce  dimensionless  energies  by  z,  = 
e«'/l  P  |1/2.  x!  =  «//|  P  |1/2.  xj  =  ?//|  0  |!/2,  z0  =  hw0/|  P  |1/2-  The  final  state 
Xf  can  be  obtained  by  applying  the  direct  technique  [2.31]  to  Eq.  (2.10).  Once 
*i,  XO)  xj  and  0  are  given,  x/  is  then  found  by  generating  random  numbers  r, 
evenly  distributed  between  (0,1),  and  using: 

a  Q  Y 

XJ  =  (x,  -  x0  -  5T|-^-|)  -  tg[aretg[xi  -  x0  -  *  j-^-j)  -  -j .  (2.13) 

Figures  2.1  -  2.3  show  the  joint  spectral  densities  as  obtained  from  Eq.  (2.7) 
(continuous  curves)  and  the  approximate  expression  (2.10)  (dashed  curves).  As 
is  seen  from  Fig.  2.1  and  2.2,  the  presence  of  ICFE  results  in  a  broadening  and 
skewing  of  the  original  delta  function.  The  skewness  depends  on  the  direction 
of  q  with  respect  to  E  and  the  broadening  can  be  so  large  that  the  final  kinetic 
energy  of  a  carrier  may  be  greater  than  the  initial  one,  even  if  a  phonon  emission 
process  has  occurred.  (Note  that  for  illustrative  purposes  in  the  figures  we  have 
chosen  the  maximum  q  available  in  the  process.)  By  decreasing  the  electric 
field  strength,  the  skewness  and  broadening  reduce  and  the  oscillations  become 
faster.  Both  Eqs.  (2.7)  and  (2.11)  recover  a  delta  function  shape  for  E  =  0. 
We  point  out  that  by  generating  104  random  numbers  according  to  Eq.(2.13) 
the  distribution  of  x/  is  found  to  coincide  with  Eq.  (2.11). 

To  introduce  ICFE  into  a  Monte  Carlo  program  we  suggest  the  following 
procedure: 

If  k,  is  the  wave  vector  before  scattering  then: 

(a)  generate  the  direction  of  the  final  k /  assuming  isotropic  distribution; 

(b)  evaluate  the  magnitude  of  kj  as  in  the  absence  of  ICFE; 

(c)  evaluate  0  using  the  kf  as  in  the  absence  of  ICFE; 

(d)  generate  xj  from  Eq.  (2.13)  and  determine  k/  in  the  presence  of  ICFE 
accordingly. 
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The  above  procedure  contains  two  weak  points.  First,  the  kj  used  to 
calculate  (3  is  not  consistent  with  the  final  k /  and  second,  we  must  fix  Xf 
appropriately.  The  first  point  leads  to  an  incorrect  determination  of  the  shift 
of  due  to  the  electric  field.  However,  since  this  shift  is  proportional  to 

q the  error  thus  introduced  should  be  practically  irrelevant. 

The  second  point  is  more  serious.  Allowing  S/  to  approach  infinity  leads 
to  carrier  run-away  and  thus  to  a  diverging  mean  kinetic  energy,  as  evinced  by 
our  simulations.  Since  this  point  is  common  to  CB  as  well,  we  will  comment  on 
it  below. 

2.1. 4.2  Collisional  Broadening 

We  shall  confine  our  interest  to  the  case  where  the  scattering  rate  is  a  function 
of  hu  only,  because  this  is  the  only  case  that  can  be  handled  analytically  (scat¬ 
tering  with  non-polar  optical,  intervalley,  and  acoustic  phonons  under  elastic 
and  energy-equipartition  approximations  are  the  corresponding  cases  of  inter¬ 
est  [2.32]).  By  taking  the  self  energy  in  the  lowest  order  in  the  electron-phonon 
coupling  [2.17]  we  find  [2.30]: 

(tfTf)  (Jly  o 

(y  -  Vo)1/2(y  -  2 y0)1/2%i  -  y0)fl(y/  -  y0) 

[(y  -  y»)2  +  (y  -  yo)][(y  -  yo  -  y/)2  +  (y  -  2y0)] 

|  ^y/0(y»  -  yo)fl(yo  -  y/)fl(y/) 

(yo  +  y/  -  yi)1  +  yj 

+  +  y /  -  y.-)0(yo  -  y.)%o  -  y /) }  (2.14) 

where  wc  use  dimensionless  energies,  y  =  hw/ 7*,  =  e./'y2,  y/  =  C//72, 

yo  =  ftwo/72,  hw 0  is  the  optical  phonon  energy.  The  parameter  7  represents 
the  coupling  strength  in  energy  units.  It  is  related  to  the  usual  deformation 


potential  constant  DtK  [2.32]  by  [7  =  m3^2(DtK)2 /(25/57rpft2wo)],  where  p  is 
the  density  of  the  material. 

The  joint  spectral  density  of  Eq.  (2.14)  is  positive  definite  but  it  exhibits  a 
long  tail  for  positive  energies,  where  it  decays  as  fy  3^2.  An  analytical  expression 
which  approximates  Eq.  (2.14)  accurately  for  ya- '»  yo  (typically  y<  >  5yo),  and 
controls  the  low  and  high  energy  tails  (K(et,e/)  =  0  for  yy  <  0  and  yy  >  t/y, 
V/  =  tfh2)  is  given  by: 


2 »; 


i/i 


with 


i2  (y.  -  if  -  yo)2  +  4y, 


„  i  .  ,yi  —  yo,  .  ,y.  -y/-yo,._i 
B  =  [«rety(— JT^-)  -  arcig(-  ^  )] 


(2.15) 


(2.16) 


We  find  the  final  state  y /,  once  y,-,  yo  and  yy  are  given,  by  generating  yy  as: 


y/  =  y.  -  yo  +  2y,1/a  +  (£  -  aretg( (2.17) 

where  r  is  a  random  number  uniformly  distributed  between  (0,1). 

Figures  2.4  -  6  show  the  joint  spectral  densities  as  obtained  from  Eq.  (2.14) 
(continuous  curves)  and  the  approximate  expression  (2.15)  (dashed  curves).  As 
a  general  trend  the  approximation  improves  for  increasing  values  of  the  initial 
kinetic  energy  (see  Figs.  2.4  and  2.5).  As  can  be  seen  from  Figs.  2.5  and 
2.6,  the  broadening  increases  for  increasing  7  values,  and  it  gives  rise  to  carrier 
run-away  as  in  the  case  of  ICFE.  We  point  out  that  by  generating  104  random 
numbers  according  to  Eq.  (2.17)  the  distribution  of  yy  is  found  to  coincide  with 
Eq.  (2.16). 

To  introduce  CB  into  a  Monte  Carlo  program  we  suggest  the  following 
procedure. 

If  k,-  is  the  wave  vector  before  scattering,  then: 
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(i)  generate  the  direction  of  the  final  k /  assuming  isotropic  distribution; 

(ii)  generate  y /  from  Eq.  (2.17)  and  determine  kj  accordingly. 

The  above  procedure  requires  an  appropriate  choice  for  y j  just  as  was  the 
case  of  ICFE. 

2.1.5  Numerical  Calculations 

To  investigate  the  main  features  of  ICFE  and  CB  we  have  performed  Monte 
Carlo  simulations  using  104  electrons  for  the  simple  model  semiconductor  char¬ 
acterized  by  the  three  parameters  m  =  0.3mo  (mo  is  the  free  electron  mass), 
hwo  ~  40meK  and  -y2  =  1.1  meV.  The  choice  of  these  values  can  be  considered 
as  typical  for  several  cubic  semiconductors.  As  a  plausible  choice,  we  have  taken 
the  cut-off  energy  in  dimensionless  units  as  x/  =  2x,  for  ICFE  and  y j  =  2y,  for 
CB,  respectively. 

Figures  2.7  to  2.10  display  the  energy  distribution  functions  as  obtained 
in  the  presence  of  ICFE  alone  and  CB  alone  for  various  different  electric  field 
strengths.  For  the  sake  of  comparison,  we  have  also  given  the  results  for  the 
semiclassical  case  (SC)  obtained  by  using  a  delta  function  for  the  joint  spectral 
density.  (In  order  to  facilitate  a  detailed  comparison  all  the  curves  show  the 
direct  Monte  Carlo  results). 

For  the  high  electric  fields  (500  and  100  kV/cm)  the  presence  of  ICFE  and 
CB  is  found  to  strongly  increase  the  number  of  carriers  in  the  high  energy  tail 
of  the  electron  distribution  with  respect  to  the  SC  case.  We  have  also  found 
that,  in  general,  CB  is  responsible  for  an  increase  in  the  population  of  the  lower 
energy  region  of  the  distribution  function,  while  ICFE  evidences  an  opposite 
effect.  At  these  high  fields,  the  SC  model  is  characterized  by  a  heated  Maxwell- 
Boltzmann  distribution,  as  expected  within  the  quasi  elastic  regime  [2.32],  In 
particular,  we  notice  that,  for  the  present  choice  of  parameters,  ICFE  gives  a 
slightly  larger  effect  than  CB. 
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For  the  low  electric  fields  (10  and  5  kV/cm),  the  calculations  show  that  the 
effect  of  ICFE  and  CB  becomes  negligible  (see  Figs.  2.9  and  2.10),  apart  for 
some  minor  effects  of  ICFE  in  the  lower  energy  tail  of  the  distribution  probably 
due  to  the  approximate  expression  (2.13)  used.  At  these  low  fields,  the  SC 
model  approaches  the  streaming  motion  regime’  [2.32],  where  the  distribution 
function  sharply  drops  to  zero  for  kinetic  energies  larger  than  the  optical  phonon 
energy. 

Figure  2.11  shows  the  distribution  of  the  energy  difference  between  states 
before  and  after  a  scattering  event  ( Ae  =  e,  —  e /).  With  respect  to  the  SC  case, 
for  which  Ae  =  ftw o,  both  ICFE  and  CB  exhibit  a  broadened  distribution  which 
reflect  the  joint-spectral-density  model  used  in  the  simulation.  In  particular, 
the  two  peaks  for  the  case  of  ICFE  are  correlated  to  the  skewness  associated 
with  the  forward  and  backward  directions  of  q  with  respect  to  E. 

Table  2.1  summarizes  the  values  of  the  mean  kinetic  energy,  the  drift  veloc¬ 
ity,  and  the  maximum  kinetic  energy  achievable  by  a  carreir  during  a  simulation 
for  the  different  cases.  At  low  fields  SC  and  CB  give  the  same  results,  within  the 
numerical  uncertainty,  while  ICFE  gives  a  value  which  deviates  for  about  20%, 
with  a  larger  mean  energy  and  a  smaller  drift  velocity  than  the  SC  case.  These 
deviations  should  be  attributed  to  the  scarce  reliability  of  the  model  spectral 
density  at  the  lowest  energies.  At  high  fields  the  mean  energy  for  ICFE  and  CB 
is  significantly  larger  (over  50%  at  500  kV/cm)  than  for  the  SC  case,  reflecting 
a  higher  population  of  the  high  energy  tail  of  the  distribution.  Conversely,  the 
drift  velocity  change  is  limited  to  slight  deviations,  about  15%.  However,  an 
increase  in  the  value  of  the  high  energy  cut-off  ?/  of  the  joint  spectral  density 
has  been  found  to  be  responsible  for  a  systematic  increase  in  the  mean  kinetic 
energy.  As  a  consequence,  the  maximum  kinetic  energy  achievable  by  a  carrier 
during  a  simulation  also  increases  together  with  its  numerical  uncertainty,  and  it 
is  clear  that  the  computer  time  necessary  for  a  simulation  becomes  unacceptably 
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long.  This  effect  is  shown  in  Table  2.2  for  the  case  of  ICFE  at  E  =  100 kV/cm. 
Analogous  results  have  been  found  for  the  case  of  CB  as  well.  Clearly  this  in¬ 
crease  of  the  mean  kinetic  energy  is  associated  with  a  larger  population  of  the 
high  energy  tail  of  the  distribution  function.  Therefore,  we  conclude  that  the 
choice  of  the  high  energy  cut-off  is  important  for  a  quantitative  evaluation  of  the 
effects  of  both  ICFE  and  CB.  This  is  a  crucial  point  for  the  present  modeling 
and,  in  our  opinion,  it  has  affected  all  previous  calculations  [2.16  -  24].  Thus, 
physical  justifications,  which  are  outside  the  scope  of  this  paper,  are  needed  to 
overcome  this  drawback. 

2.1.6  Conclusions 

We  have  presented  a  theoretical  framework  which  enables  one  to  account 
for  genuine  quantum  effects  for  the  case  of  hot-electron  transport  in  semiconduc¬ 
tors.  The  main  quantities  of  interest  are  the  spectral  density  and  the  associated 
joint  spectral  density  which  can  include  ICFE  as  well  as  CB.  After  providing 
the  general  theory  ,  analytical  models  for  these  quantities  have  been  obtained 
in  some  cases  of  interest.  When  quantum  effects  are  neglected,  that  is  carriers 
behave  as  a  free  particle  between  two  scatterings,  the  semiclassical  Boltzmann 
picture  with  the  collision  term  determined  from  the  golden  rule  is  recovered. 
We  have  presented  two  algorithms,  compatible  to  standard  Monte  Carlo  codes, 
for  estimating  quantum  effects,  such  as  intra-collisional  field  effects  and  colli- 
sional  broadening,  in  semiconductor  high-field  transport.  The  basic  feature  of 
these  algorithms  is  to  model  joint  spectral  densities  replacing  the  Dirac  function 
which,  within  a  semiclassical  approach,  conserves  the  kinetic  energy  within  a 
scattering  event.  The  methodology  chosen  by  us  has  been  given  with  details,  so 
that  it  may  constitute  a  first  standard  which  can  be  easily  reproduced  and/or 
improved  by  other  research  groups.  The  broadening  of  the  final  energy  after 
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a  scattering  event,  which  is  associated  to  these  models,  yields  in  general  a  sig¬ 
nificant  increase  of  the  population  in  the  high  energy  region  of  the  distribution 
function.  Accordingly,  the  concept  of  a  broadening  assisted  run-away  effecty 
has  been  introduced.  However,  because  of  their  long  tails  at  high  energies, 
these  joint-spectral  density  models  need  an  ”ad  hoc"  truncation  to  avoid  an 
unpractical  (and  unphysical)  blowing  up  of  the  mean  kinetic  energy.  Because  of 
this,  in  turn,  the  quantitative  results  may  depend  on  the  choice  of  the  empirical 
cut-off  energy.  To  overcome  this  drawback,  which  has  received  some  attention 
in  the  recent  literature  [2.33],  more  appropriate  physical  models,  such  as  relax¬ 
ing  the  simple  effective  mass  model  for  the  density  of  states  [2.34],  should  be 
introduced.  We  shall  address  this  issue  in  our  further  research. 

2.2  A  Monte  Carlo  Approach  to  the  Solution  of  the  Liouville-von 

Neumann  equation  for  quantum  transport 

2.2.1  Introduction 

In  the  last  decades  the  Monte  Carlo  method  [2.31]  has  proved  to  be  a 
formidable  tool  for  the  solution  of  charge-transport  problems  in  semiconductors. 
It  provides  a  simple  way  to  solve  the  semiclassical  Boltzmann  transport  equation 
once  the  band  structure  and  the  scattering  mechanisms  are  sufficiently  well 
known.  From  its  original  formulation  [2.35]  it  has  been  greatly  improved  in 
efficiency  and  it  has  been  extended  to  cover  a  large  variety  of  problems  such  as 
space  or  time-dependent  systems,  interacting  electrons,  degenerate  statistics  or 
even  device  modelling. 

However,  the  very  fast  evolution  of  miniaturization  techniques  for  semicon¬ 
ductor  technology  is  leading  very  rapidly  towards  experimental  conditions  where 
typical  lenghts  are  of  the  order  of  the  De  Broglie  wavelenght  of  the  carriers,  and 
typical  times  are  comparable  with  carrier  relaxation  times  [2.11.  In  particular, 
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in  modern  laser  spectroscopy  a  time  resolution  has  been  achieved  of  the  order 
of  101  femtoseconds.  If  the  system  is  observed  at  a  time  of  this  order  after  it 
has  been  "prepared”,  transitions  may  be  observed  that  would  not  be  allowed 
by  energy  conservation.  In  other  words,  the  quantum  interference  phenomena 
that  produce  the  S  of  energy  conservation  are  not  yet  complete  at  so  short 
observation  times. 

It  is  clear  that  the  classical  transport  theory,  based  on  Boltzmann  equation, 
is  not  adequate  for  the  description  of  the  physical  processes  that  are  taking  place 
on  this  time  scale.  In  fact,  semiclassical  transport  is  based  on  the  hypothesis, 
among  others,  that  each  scattering  event  is  completed  when  the  next  one  starts. 
For  the  validity  of  such  an  assumption  it  is  necessary  that  the  coupling  between 
electrons  and  the  scattering  agents  is  sufficiently  weak  so  that  a  first-order 
perturbation  theory  can  be  applied,  and  this  must  be  done  in  the  limit  of 
"completed  collisions”  so  that  energy  conservation  holds  at  each  interaction 
process. 

Several  possible  approaches  have  been  presented  for  attaching  the  problem 
of  quantum  charge  transport  [2.1]. 

A  new  Quantum  Monte  Carlo  (QMC)  procedure  has  been  recently  devel¬ 
oped  [2.36]  for  the  solution  of  the  Liouville  equation  for  ihe  electronic  density 
matrix  in  semiconductors.  The  principles  of  the  method  will  be  summarized  in 
Sects.  2.2.2  and  2.2.3. 

In  Sec.2.2.4.I  we  present  an  analysis  of  the  first  two  perturbative  corrections 
in  presence  of  an  arbitrary  high  electric  field,  and  in  Sec.  2.2.4. 2  the  method  is 
applied  to  the  problem  of  energy  relaxation  of  photoexcited  electrons  in  GaAs 
in  absence  of  electric  field. 
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In  the  semiclassical  limit  our  quantum  equation  recovers  the  semiclassieal 
Boltzmann  equation,  and  the  numerical  procedure  results  to  be  a  new  formu¬ 
lation  of  the  traditional  Monte  Carlo  procedure  [2.37].  The  semiclassical  limit 
of  the  procedure  is  described  in  Sec.  2.2.5. 

The  method  allows  to  evaluate  the  electronic  density  matrix  as  a  func¬ 
tion  of  time  without  any  assumptions  on  the  intensity  and  the  duration  of  the 
electron-phonon  interaction.  The  quantum  equation  is  solved  through  a  ran¬ 
dom  generation  of  all  possible  quantum  interactions  at  the  various  perturbative 
orders,  in  the  same  way  as  the  usual  Classical  Monte  Carlo  (CMC)  generates 
classical  scattering  events  [2.36]. 

2.2.2  Physical  System  and  Theoretical  Approach 

In  order  to  study  the  properties  of  charge  transport  in  a  quantum  scheme, 
let  us  consider  an  ensemble  of  electrons  in  a  semiconductor  crystal,  coupled  to 
the  phonon  gas.  Carriers  are  assumed  to  be  not  interacting  with  each  other,  so 
that  the  interaction  of  one  carrier  with  the  phonons  will  represent  the  behaviour 
of  the  whole  electron  gas.  The  electron  band  structure  is  introduced  in  the 
effective-mass  approximation,  with  a  simple  spherical  and  parabolic  band. 

The  Hamiltonian  of  the  system  is  given  by 

H  =  Ht  +  He  +  Hp  +  Hep  (2.18) 

where  H,  —  -  V3  is  the  term  corresponding  to  an  electron  in  a  perfect  crystal 
(m~  effective  mass);  He  =  eE  ■  r  describes  the  electric  field  turned  on  at  t=0; 
Hp  =  YL  h0Jqalaq  describes  the  free-phonon  system  in  the  second-quantization 
formalism  (oj  and  a,  axe  the  creation  and  annihilation  operators  of  a  phonon 
mode  q).  The  electron-phonon  interaction  Hamiltonian  Htp,  turned  on  at  t  =  0, 
depends  on  the  scattering  mechanisms  included  in  the  model. 
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We  have  not  explicitely  introduced  any  interactions  among  phonons  and 
between  phonons  and  the  thermal  bath.  In  the  numerical  procedure,  however, 
we  will  assume  that  these  interactions  can  mantain  an  equilibrium  phonon  pop¬ 
ulation  during  the  evolution  of  the  system. 

In  order  to  work  in  the  interaction  representation,  we  use  the  set  of  basis 
functions 

I  k,  {*»,},  t)  =  -Lt'W)0.(-7:^(Mt)).$({Ri})t)  (2.19) 

that  include  the  time  dependence  due  to  the  unperturbed  Hamiltonian  Hc  + 
Hp  +  He-  They  are  direct  products  of  electronic  accelerated  plane  waves,  with 
k(t)  =  k0  —  *^t,  normalized  to  1  over  the  crystal,  and  the  phonon  states 
$({n?})  with  nq  phonons  in  mode  q  with  frequency  uiq.  The  state  |  'll)  of  the 
system  can  be  expanded  over  this  set  as: 

I  *>  =  £  £  e(Mn,},0  |k,{n,},t>  (2.20) 

k  <«,) 

If  we  now  consider  to  the  density  matrix  of  the  system  in  the  representation 
of  the  set  in  Eq.(2.19): 

*>(k,{n,},k',{n^},t)  =  (e(k,  {n,},t)e*(k',  (n'},t)),  (2.21) 

the  Liouville-Von  Neumann  equation  that  describes  its  time  evolution  contains 
only  the  perturbation  Hamiltonian: 

ih~p(X,X',t)  =  [He  p,p\(X,X'M),  (2.22) 

where  we  have  used  the  symbolic  compact  notation  X  =  (k,  {nq}) 

A  formal  integration  leads  to 

p(X,X',t)  =/>(*,  X',0)+  f*  dhlH^pKX.X^h),  (2.23) 

Jo 
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where  Hep  =  ^Hcf-  If  we  are  interested  in  the  evaluation  of  expectation  values 
of  electron  quantities  which  are  diagonal  in  the  electronic  part  of  the  states  in 
Eq.(2.19),  we  can  focus  our  attention  on  the  diagonal  elements  p(X,  t)  of  p. 

Furthermore,  we  will  assume  a  diagonal  initial  condition  for  p  decoupled 
in  electron  and  phonon  coordinates.  This  is  justified  by  the  fact  that  the  inter¬ 
action  is  turned  on  at  t=0.  The  electronic  part  is  taken  as  some  distribution 
function  /o,  while  the  phonon  part  is  assumed  as  the  probability  Pcq{{nq})  of 
finding  each  mode  q  occupied  by  nq  phonons  at  equilibrium: 

p(kf{n,},0)  =  /„(k)-/V,({n,})  (2-24) 

A  perturbative  expansion  of  Eq.(2.23)  is  easily  obtained  by  iterative  sub¬ 
stitution  of  its  right-hand  side  into  the  equation  itself: 

p(x,t)  =  p(x,o)  +  /*<#i(M*i).p(o)IPQ+ 

Jo 

+  fdh  o)H(*)  +  ...= 

Jo  Jo 

=  p{o){x,t)  +  ap(1)(a:,o  +  a  pw(x,t)  + ...  (2.25) 

The  zero-order  term  in  the  expansion  corresponds  to  the  case  of  no  coupling 
between  electrons  and  phonons,  and  it  is  equal  to  the  initial  condition  for  p. 

In  order  to  proceed  with  the  theory  we  note  that  the  perturbation  Hamil¬ 
tonian  has  the  form 

Hep  =  £/•(«!)<-,«*-  -  r}  =  B.t  +  Hem,  (2.26) 

V 

where  Hab  and  Htm  refer  to  phonon  absorption  and  emission,  respectively.  The 
matrix  element  of  Hrp  between  X  and  X‘  contains  only  the  mode  q  related  to 
k  and  k'  by  momentum  conservation  and  it  is  different  from  zero  only  if  the 
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number  of  phonons  in  the  mode  q  is  changed  by  a  unity  going  from  X  to  X', 
since  it  contains  only  linear  terms  in  aq  e  aj. 

The  explicit  form  of  the  first-order  correction  is: 

Ap(1)(x,o  =  f*dtl{)imp[x,x,,tl)p^[x,to)  -  jfv(xtx,f*1)p<0>(*,o)}. 

Jo 

(2.27) 

From  what  we  have  seen  about  the  matrix  elements  of  Hcp,  it  is  clear  that  the 
above  first-order  term  gives  no  contribution  to  the  diagonal  elements  of  p  since 
we  assumed  diagonal  initial  conditions.  It  can  be  shown  that  the  same  is  true 
for  all  odd-order  corrections. 

The  second-order  term  can  be  rewritten,  using  the  property 

Xep(X,X',t)  =  -Kf(X',X,t),  (2.28) 

as 

Ap(2)(X,t)  =  Y.  !  dt'  /  '  dt2{Hep(X,X\h)Xcp{X',XM)p{X,0) 

x,  Jo  Jo 

+x;p(x1x',tl)x;p(x',x,h)P(x,o)  +  x;p(x,x',tl)xep(x,x',t2)p[x',o) 
+xtp(x,x',t1)x;p(x,x',t2)p(  a'.o)}  (2.29) 

There  is  a  simple  and  useful  way  of  reading  the  above  equation.  At  t  =  0  the 
two  arguments  of  p  are  equal;  by  application  of  Xcp  (or  Xqp)  the  first  (or  the 
second)  argument  of  p  is  changed  from  the  second  argument  of  Xep  to  the  first 
one;  at  l  the  two  arguments  of  p  are  again  equal  (to  X).  Since  each  application 
of  Xcp  changes  the  phonon  state  of  one  unity,  in  order  to  start  from  a  diagonal 
element  and  end  up  to  another  diagonal  element  a  mode  q  absorbed  (or  emitted) 
by  one  argument  is  to  be  absorbed  (or  emitted)  also  by  the  other  argument  or 
reemitted  (or  reabsorbed)  by  the  same  argument.  Using  the  language  of  the 
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field  theory,  we  refer  to  the  first  kind  of  processes  as  to  "real"  emissions  and 
absorptions,  while  the  other  ones  are  "virtual  processes”. 

Thus,  the  contributions  to  be  included  in  the  second  order  in  Utp  are  those 
illustrated  in  Fig.  2.12. 

With  the  above  interpretation  it  is  very  simple  to  generalize  the  results 
to  higher-order  terms  of  the  perturbative  expansion.  For  example,  Fig. 2. 13  is 
the  diagramatic  representation  of  the  following  contribution  to  the  fourth-order 
term: 

Km  {X,  X\t1))ial>{X,  X"\  t3)Xm(Xm,  X",t3)){tm(X",  X\  t4)p(X',  0)  (2.30) 

Each  "process"  in  these  graphs  corresponds  to  a  single  scattering  event  in 
classical  transport.  To  recover  the  classical  golden  rule  it  should  be  necessary 
to  integrate  one  of  the  two  times  of  the  process,  without  interfering  with  other 
processes,  over  an  interval  large  enough  to  obtain  the  6  of  energy  conservation. 

From  the  analysis  of  these  graphs  of  the  perturbation  expansion  for  p,  we 
can  focuse  our  attention  to  some  important  aspects  of  the  quantum  description. 
The  quantum  transitions  have  a  finite  duration,  during  which  carriers  experience 
the  action  of  the  electric  field  (ICFE).  This  effect  is  obviously  not  present  in 
the  semiclassical  description  where  collisions  are  point-like.  Only  after  a  certain 
time  we  obtain  again  a  diagonal  state  that  corresponds  to  a  semiclassical  state 
of  the  system;  during  the  interaction  the  system  is  in  a  quantum  state  given 
by  a  superposition  of  k  states  which  does  not  correspond  to  any  semiclassical 
situation.  Furthermore,  while  one  process  is  happening,  a  new  one  can  start, 
giving  rise  to  multiple  collisions  and  vertex  corrections. 

Let  us  now  introduce  the  reduced  electronic  density  matrix,  which  is  the 
quantity  of  interest  in  charge-trasport  phenomena  and  it  can  be  considered  the 
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quantum  analogue  of  the  classical  distribution  function.  It  is  defined  as  the 
contraction  over  phonon  states  of  the  total  density  matrix: 

p(<)(k,k',f)  =  P(ki(nJ>k'. {«?}>*)  (2.31) 

The  evolution  equation  for  p W  still  contains  the  dynamical  variables  of  the 
many-body  system.  In  fact  it  is  not  possible  to  obtain  a  closed  equation  for 
pW  by  taking  the  trace  of  Eq.(2.23)  for  the  full  density  matrix  since  the  trace 
operation  does  not  commute  with  the  interaction  Hamiltonian.  The  reduction 
of  the  total  density  matrix  to  the  electron  density  matrix  can  instead  be  easily 
performed  in  the  numerical  QMC  procedure  introduced  in  the  following  section 
in  order  to  evaluate  the  series  in  Eq.(2.25). 

2.2.3  Numerical  Procedure 

The  numerical  QMC  algorithm  devised  for  the  solution  of  Eq.(2.22)  is 
essntially  based  on  random  generations  of  all  possible  processes  associated  with 
the  different  perturbative  corrections. 

From  the  fundamentals  of  the  Monte  Carlo  technique  [2.38]  it  is  known  that 
for  the  evaluation  of  a  sum  S  =  x,  ,  one  can  consider  the  estimator  ^ ,  where 
Pi ’s  are  arbitrary  probabilities  between  zero  and  one  normalized  to  unity,  and 
average  it  over  random  selections  of  the  i  index,  made  with  probability  />,  This 
procedure  is  here  used  to  obtain  an  estimate  of  Eq.(2.25):  random  selections 
with  suitable  probabilities  determine 

i.  the  order  of  the  perturbative  correction  to  be  estimated; 

ii.  one  of  the  possible  contributions  to  the  corresponding  integrand; 

iii.  the  wavevectors  q  of  the  phonons  involved  in  the  quantum  interactions 
of  that  contribution.  Due  to  momentum  conservation  of  the  Hep  matrix 
elements,  these  selections  determine  the  argument  kj„  of  p  at  t=0. 


The  quantity: 


dtnHep{ti)...X'P{tn) 


—  0)t 


where  P  is  the  total  probability  of  all  the  selections  that  have  been  made,  is 
then  averaged  over  many  generations  and  it  gives  an  estimate  of  p  at  time  t. 

In  this  way  we  obtain  the  density  matrix  and  the  average  values  of  the 
physical  quantities  of  interest  for  our  system  through  a  random  generation  of 
quantum  processes  as  we  obtain  the  carrier  distribution  function  and  transport 
quantities  from  a  random  choice  of  carrier  histories  in  the  traditional  Monte 
Carlo  technique. 

There  is  however  a  crucial  point  that  requires  to  be  analised,  which  is 
how  we  perform  the  average  over  the  phonon  variables.  As  we  already  pointed 
out,  we  assume  an  initial  condition  for  p  which  is  the  product  of  an  electron 
distribution  function  times  the  equilibrium  phonon  distribution;  furthermore 
the  interaction  Hamiltonian  is  linear  in  the  creation  and  annihilation  operators 
of  the  modes  q. 

In  the  numerical  procedure  we  saw  that  a  sequence  of  process  is  gener¬ 
ated  starting  from  an  initial  state  (k,„,  which  terminates  on  the  state 

(k,{n,},t).  This  sequence  corresponds  to  a  sequence  of  phonon  wavevectors  q 
which  are  absorbed  or  emitted  on  the  first  or  on  the  second  index  of  the  density 
matrix.  We  can  assume  that  each  phonon  mode  is  chosen  only  once  in  a  given 
process,  i.e.  the  electron  always  interacts  with  an  equilibrium  phonon  bath.  If  a 
phonon  q  is  absorbed  in  a  certain  transition  and  the  corresponding  occupation 
number  at  time  t  is  nq,  the  occupation  number  of  the  initial  state  is  nq  +  1,  and 
the  interaction  hamiltomnian  l/ep  must  contain  the  factor  \/ (riq  +  1)  because  it 
contains  the  operator  aq.  In  order  to  finish  the  process  Pep  must  act  two  times, 
and  the  numerical  estimator  will  contain  the  multiplicative  factor  nq  +  1. 
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In  the  same  way,  if  we  have  the  emission  of  a  mode  q,  and  the  corresponding 
final  occupation  number  is  nq,  then  the  initial  occupation  of  this  mode  must  be 
nq  —  1,  and  the  numerical  estimator  will  contain  a  factor  nq. 

At  the  end  of  the  random  generation  we  obtain  for  the  2n-order  correction 
Ap('K2n)(k,t)  an  estimator  of  the  form 

Ap<‘K2">(k,t)  =  (constant)  ■  fa( kt„)  •  £ {[[Nq.  •  f[^K  J)  (2-33) 

{"«}  «'  « 

where  n1in  is  the  occupation  number  of  the  mode  q  in  the  initial  state  and 
Nq-  is  the  multiplicative  factor  generated  by  the  interaction  hamiltonian  for 
the  phonon  modes  q*  involved  into  the  generated  sequence.  The  average  (...)  is 
performed  over  several  generations  of  terms  for  a  given  perturbative  order  2n. 

If  a  phonon  q  is  not  chosen,  then  the  sum  over  all  possible  occupation 
numbers  nq  of  P(nq)  must  be  equal  to  unity: 

£PKn)  =  l  (2.34) 

»« 

and  this  factor  does  not  contribute  to  the  product  in  Eq.  (2.33). 

If  instead  we  consider  a  q*  chosen  in  the  random  generation,  we  have  two 
possible  case.  In  the  case  of  absorption: 

£<K*  +  i)^K*  + 1))  —  ^*)bok  (2.35) 

V 

where  nq~)Bote  is  the  equilibrium  occupation  number  given  by  the  Bose  distri¬ 
bution. 

In  the  case  of  emission: 

^((n,.)P(n,<  -  1))  =  nq.)Bo.e  +  1  (2.36) 
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In  both  cases  we  have  finally  to  use  the  Bose  distribution  for  the  evaluation 
of  the  term  without  introducing  any  approximations  on  the  electron-phonon 
coupling,  but  for  the  assumption  that  the  phonon  gas  is  constantly  kept  in 
equilibrium  conditions. 

2.2.4  Results 

2.2.4. 1  Analysis  of  First  and  Second  Orders  in  Presence  of  Electric  Fields 

The  second  order  perturbative  correction  in  the  interaction  hamiltonian 
involves  only  one  of  the  processes  discussed  above.  The  explicit  form  of  the 
corresponding  integrals  in  Eq.(2.25)  is  of  the  type: 

f  dti  f 

Jo  Jo 

a=_L_(k,  — kll.E 

*“;("/-■» i-"i>  P-37) 

Here  i  and  f  refer  to  the  interacting  initial  and  final  electron  states  and  uiq  is  the 
phonon  frequency.  This  simple  expression  allows  a  direct  analytical  integration 
in  terms  of  Fresnel  integrals.  The  final  expression  of  the  second  order  correction 
after  integration  is: 

JL{[C(^±1)  -  C(  V  +  (S(*±*)  -  S(£)]2}  (2.38) 

In  the  above  equations  the  coefficient  a  accounts  for  the  effect  of  the  field  dur¬ 
ing  the  finite  duration  of  the  collision,  while  b  is  related  to  the  energy  of  the 
quantum  states  involved  in  the  transition. 

This  formulation  allows  us  to  investigate  separately  the  new  features  in¬ 
troduced  by  the  quantum  treatment  of  the  interactions.  In  particular,  we  can 


turn  off  the  ICFE  by  neglecting  the  effect  of  the  field  between  two  vertices  of 
the  process. 

Results  for  the  first  two  perturbative  terms  have  been  obtained  starting 
from  equilibrium  conditions  for  the  electron  and  phonon  system,  and  both 
electric  field  and  electron-phonon  coupling  are  turned  on  at  t=0.  In  Fig.2.14 
we  compare  the  absolute  value  of  the  second  order  correction  as  given  from 
Eq.(2.38)  with  the  same  contribution  without  ICFE.  In  the  same  figure  we 
report  also  the  corresponding  classical  contribution  to  the  perturbative  expan¬ 
sion  of  the  classical  integral  equation  for  the  distribution  function  obtained 
from  the  Boltzmann  equation  [2.37].  This  classical  correction  comes  from  "one- 
scattering”  trajectories  and  the  corresponding  contribution  is  always  negative, 
due  to  the  prevailing  scattering  out.  These  numerical  results  have  been  obtained 
for  a  simple-model  semiconductor  (relative  effective  mass  m*=. 295,  crystal  den¬ 
sity  d=2.‘33  2/cm3,  one  optical  phonon  scattering  with  Tpk  —  350 K,  and  cou¬ 
pling  constant  D=2.5109eVcm~1  ).  The  working  conditions  (T=20  K,  E=150 
kV/cm,  t=5.xl0-1'*  s)  have  been  chosen  in  such  a  way  that  quantum  effects 
can  be  easily  detectabie. 

When  we  neglect  the  ICFE,  we  get  a  higher  effect  of  the  "quantum  one- 
collision”  trajectories.  In  fact  the  ICFE,  by  changing  the  energy  of  the  carrier 
during  the  collision,  reduces  the  efficiency  of  the  scattering  since  it  reduces  the 
time  of  positive  interference  which  occurs  when  the  energy  difference  between 
initial  and  final  states  is  equal  to  the  phonon  energy  [2.39].  The  classical  con¬ 
tribution  is  even  higher;  such  effect  can  be  interpreted  by  noting  that  the  field 
is  so  high  that  even  at  such  short  times  the  electrons  can  reach  enough  energy 
for  classical  instantaneous  phonon  emissions,  while  the  time  is  so  short  that 
quantum  transitions,  with  finite  durations,  are  not  yet  fully  developed. 
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An  improvement  in  the  efficiency  of  the  method  can  be  obtained  also  at 
higher  orders  by  integrating  over  time  every  other  vertex  between  the  two  ad¬ 
jacent  verteces.  The  result  is  again  expressed  in  terms  of  Fresnel  integrals. 

If  we  allow  only  processes  that  do  not  overlap  (the  two  verteces  of  one  given 
processes  correspond  to  adjacent  times)  we  neglect  multiple  collisions,  and  by 
comparison  we  may  analyse  their  effect. 

Finally  we  may  neglect  the  effect  of  the  field  in  the  interference  exponentials 
and  analise  in  this  way  the  ICFE  in  the  two-scattering  trajectories.  In  this  last 
case  the  integration  of  one  vertex  leads  to  the  functions  like  which  would 

lead  to  the  delta  of  energy  conservation  for  large  times;  in  our  case  however 
the  completion  of  the  transition  is  not  necessarily  reached  if  the  time  interval 
considered  is  very  short. 

We  show  in  Fig.2.15  the  quantum  corrections  of  the  fouth  perturbative 
order  (  which  is  now  larger  than  zero)  compared  with  the  classical  two-scattering 
contribution  for  the  same  model  as  in  Fig.2.14.  The  quantum  result  is  lower 
than  the  classical  term  as  for  the  second  order  contribution,  due  again  to  the 
short  time  available  for  the  interaction.  When  we  allow  only  separate  collisions 
(Fig.2.15c),  the  corresponding  term  appears  somewhat  higher.  This  result  can 
be  due  to  the  fact  that  the  two  separate  collisions  allow  for  a  higher  contribution 
of  transitions  which  do  not  conserve  energy,  even  though  the  phenomenon  still 
requires  further  investigation. 

If  we  turn  off  the  ICFE,  we  see  that  the  whole  curve  is  still  higher  and 
much  closer  to  the  classical  one,  as  it  happens  for  the  second  order  term. 

Finally  in  Fig.2.16  we  report  results  for  the  density  matrix  in  the  case  of 
a  more  realistic  set  of  parameters  that  gets  closer  to  a  simplified  silicon  model. 
We  changed  the  phonon  temperature  (T=450  K)  and  the  coupling  constant 
(D=.8x  109  eV/cm).  The  electric  field  is  E=15  kV/cm  and  the  time  is  t=.5  ps. 
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In  this  case  we  found  that  the  ICFE  is  much  lower  than  in  the  previous  case 
due  to  the  lower  field  strength. 

Since  the  time  is  larger,  the  effect  of  multiple  collisions  is  also  lowered,  and 
we  expect  that  this  effect  influences  higher  order  corrections. 


2.2.4.2  Quantum  Energy  Relaxation  of  Photoexcited  Carriers  in  Absence 
of  Electric  Fields 

When  electric  fields  are  absent  ( He  =0)  time  integrations  can  be  per¬ 
formed  analytically.  For  example,  the  second-order  contribution  due  to  a  real 
emission  of  a  phonon  in  mode  q  is  given  by  the  first  two  graphs  on  the  left  in 
Fig.2.12,  which  are  complex  conjugate  of  each  other.  The  integration  yelds 

2*  /‘dt,  /*'  dt2(X,h  |  Me  |  X',ti)(X,t2  |  Mb  \  X',t2y 
Jo  Jo 

=  2Z{±  fdt ,  fU  dt2  |  F(q)  |2  (n,  + 
n  Jo  Jo 


_2lF(q)|2(n,4-l) 

"■ — -d -«*(*■*)). 


(2.39) 


where  Z  indicates  the  real  part,  and  6u>  =  u){k)  —  v{k')  +  uq.  Similar  results 
are  obtained  for  other  contributions  and  for  higher-order  terms.  It  is  not  easy 
to  give  general  formula  since  the  form  of  the  result  depends  on  the  particular 
diagram  considered. 

The  method  described  in  the  previous  sections  has  been  applied  to  the  case 
of  photoexcited  electrons  in  bulk  GaAs.  The  semiconductor  model  has  been 
simplyfied  to  a  single  spherical,  parabolic  band.  The  interaction  Hamiltonian 
includes  only  polar  coupling  to  optical  phonons;  for  such  a  case 


i  m  i*.  -  £>. 


(2.40) 
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where  w0  is  the  frequency  of  the  optical  phonons,  assumed  constant,  V  is  the 
volume  of  the  crystal,  Soo  and  e0  are  the  high  frequency  and  the  static  dielectric 
constants. 

Electrons  are  generated  at  t  =  0  according  to  a  distribution  proportional 
to  exp( |  c  —  c0  |2  /(KbT{),  where  e  is  the  electron  energy  and  t0  and  Tx  are 
appropriate  constants 

Terms  up  to  fourth  order  have  been  included  in  the  numerical  computa¬ 
tions.  For  such  a  reason  the  simulation  time  has  been  kept  <  lOOps.  For  longer 
times,  higher  perturbative  orders  would  have  been  necessary. 

For  comparison  results  have  been  obtained  for  the  same  model  at  the  same 
times  in  semiclassical  transport,  using  an  Ensemble  Monte  Carlo  (EMC)  tech¬ 
nique. 

The  following  parameters  have  been  used:  m  =  0.063mo;  huajKB  =  410 K\ 
T  =  lOJff;  Ti  =  10 K;  Eoo  =  10.92;  e„  =  12.9;  ta/KB  =  1000 K 

Fig.2.17  shows  the  results  obtained  with  EMC  at  t  =  O.lps  after  excitation. 
The  echos  are  clearly  seen  corresponding  to  electrons  having  emitted  one  or  two 
optical  phonons.  About  30  %  of  the  particles  have  left  the  original  peak. 

Fig.2.18  shows  the  corresponding  result  obtained  with  quantum  transport 
theory  (note  the  scale  change).  The  initial  distribution  is  diminished  of  a  quan¬ 
tity  very  similar  to  that  of  the  classical  case.  However,  electrons  can  be  found, 
at  t  =  O.lps,  in  a  very  wide  range  of  energies,  since  energy  needs  not  to  be 
conserved.  The  secondary  peaks  are  not  yet  well  formed.  High  energies  are  in 
fact  favored  in  the  distribution  at  time  t  by  the  larger  density  of  states  and  by 
the  q~ 2  factor  in  Eq.(2.40).  Thus  the  energy  relaxation  predicted  by  quantum 
transport  theory  is  less  than  that  predicted  by  classical  transport  for  t  <  0.1  ps 
as  shown  in  Fig.2.19. 
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2.2.5  Semiclassical  Limit:  Backward  Monte  Carlo  Procedure 

The  semiclassical  limit  of  the  theoretical  approach  described  in  the  pre¬ 
vious  sections  is  obtained  when:  i)  v/e  neglect  ICFE;  ii)we  assume  that  the 
time  between  two  collisions  is  much  longer  than  the  collision  duration,  so  that 
during  each  collision  energy  is  conserved;  iii)  the  phonon  population  is  always 
mantained  at  equilibrium,  and  the  average  occupation  number  is  given  by  the 
Bose  distribution. 

Under  these  conditions  the  integral  perturbative  equation  for  the  density 
matrix  reduces  to  the  Boltzmann  equation  written  in  an  analogous  integral  form 
with  the  collision  term  expanded  at  the  same  perturbative  order. 

The  semiciassical  limit  of  the  Quantum  Monte  Carlo  gives  a  basically  new 
M.C.  method  for  the  solution  of  the  Boltzmann  equation,  which  results  as  the 
classical  limit  of  the  method  developed  for  the  solution  of  the  Liouville  equation 
for  quantum  transport  [2.36]. 

The  technique  described  here  differs  from  the  traditional  M.C.  method  in  two 
major  respects:  1).  The  occurrence  of  particular  electron  histories  with  given 
scattering  events  is  arbitrarely  selected  in  the  procedure  and  appropriately 
weighted  in  the  estimator.  2).  The  electron  state  k  at  which  the  distribu¬ 
tion  function  is  evaluated  at  time  t  is  chosen  at  the  beginning  of  the  procedure 
and  the  electron  paths  are  generated  backward  in  time  from  t  to  the  time  t=0 
of  the  (known)initial  condition.  This  second  feature  seems  to  be  not  inherent 
to  the  method  which  should  be  suitable  also  for  a  normal  forward  simulation. 
However  the  fact  that  the  value  of  k  at  which  f  is  evaluated  is  fixed  arbitrarely 
makes  the  method  particularly  appealing  for  problems  where  rare  regions  of  / 
are  of  particular  interest.  The  two  features  indicated  above  may  suggest  the 
names  "weighted  M.C.”  or  "backward  M.C.”  for  the  procedure  introduced  here. 
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2.2.5. 1  The  Method 

Let  us  start  from  the  standard  form  of  the  B.E.  For  semplicity  we  shall  con¬ 
sider  here  a  homogeneus  system,  with  a  homogeneus  and  constant  applied  field 
E,  but  the  method  should  be  easely  generalizable  to  space  and  time  dependent 
phenomena. 

|/(r,k,t)  +  /(r,k,t)  =  Jdk^kLkMr.kni)- 

J  dkiP o(k,ki)/(r,k,t)  (2.41) 

Here  P<  and  P0  indicate  the  same  transition  probability  from  a  state  k  to 
state  kj;  the  suffixes  i  and  o  are  intended  to  spcify  IN  or  OUT  scattering, 
respectively,  for  future  reference.  If  the  path  variables  k*  =  k  —  -£-t  are  used 
and  a  formal  integration  is  performed,  the  Boltman  Equation  takes  the  form  of 
an  integral  equation  that  can  be  iteratively  expanded  in  powers  of  the  scattering 
probability. 

/(k,<)  =  /0  +  /*+/J+...  (2-42) 

/°(k,t)  =  /(k(0),0) 

/‘(k,t)  =  I  j  rfk1/>  (ki,k(*1))/(ki(0),0)  — 

dklFo(k(tl),k1)/(k(0),0) 
f2(k,t)  =  Jdt  j  dkiPiOuMtO) 

J  ^*aP<(ka,k,(t8))/(ka(0),0) 

/ rf<1(2^  / / ‘*k*ft(ki(*a),k>)/(ki(0),0) 

Jdi J  dkiP0(k(<i),ki)  Jdh~ 3  J  rfk2ft(k„k(t,))/(k2(0),0) 


31 


+  /  dtl(2 J  J  ^2/>O(k(t2),i2)/(k(0),0) 

(2.43) 

Here  k (t)  is  the  value  of  k  moved  back  to  time  t  . 

The-zero  order  term,  /°,  which  corrisponds  to  the  absence  of  interaction,  is 
the  ballistic  translation  of  /  during  the  interval  (0,t),  due  to  the  action  of  the 
external  electric  field.  The  first  order  term  f1  rapresents  the  contribution  of 
electrons  that  suffered  one  scattering  event  between  t=0  and  t=t.  Here,  the  first 
integral,  or  ’IN’  term,  is  due  to  electrons  that  at  time  1 1  have  a  collision  from 
ki  to  k(ti)at  time  t.  The  second  integral  or  ’OUT’  term,  is  due  to  electrons 
that  are  scattered  out  from  the  ballistic  trajectory  and  do  not  reach  the  state 
(k,t)  and  reach  k  at  time  t.  An  interpretation  like  this  can  be  carried  on  to  all 
the  other  terms.  An  example  of  a  third  order  term  is  shown  in  Fig.  (2.20). 

Eq.  (2.20)  gives  /  formally  as  an  infinite  sum.  A  possible  M.C.  technique  for 
the  evaluation  of  a  sum  S  =  £2  X,-  is  the  following:  we  select  values  of  i  at 
random  with  arbitrary  probabilities  p,  such  that  p,  >  0  and  £)  Pi  =  1  and 
make  the  estimator  x,/pt-  so  we  have 


The  method  can  be  extended  to  the  evaluation  of  an  integral:  if  we  consider  a 
multiple  integral  like  those  in  Eq  (2.20)  then 

J  dti  J  ...  j  dtnf{tuh,...,tn)  =  v  <  />=  it"  <  /(f  >  (2.45) 

where  the  mean  value  is  evaluated  over  the  possible  choices  of  the  set 
(ti,<2,  ...,tn)  and  v  is  the  volume  of  integration. 

If  we  apply  the  above  method  for  the  evaluation  of  the  sum  in  Eq  (2.20)  the 
following  procedure  could  results:  1)  the  perturbative  order  n  of  the  term  to  be 
evaluated  is  chosen  first ;  2)  n  times,  between  t  =  0  and  the  time  t  of  simulation, 
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are  then  generated  ;  3)  then,  starting  from  the  final  state  (k,t)  we  have  to  move 
k  backward  in  time  with  free  flights  from  each  t;  to  the  previous  one  ;  4)  at  each 
*«  we  choose  between  scattering  IN  or  OUT  and  the  scattering  mechanism  (for 
example,  in  our  case  where  we  consider  only  scattering  by  one  optical  mode, 
we  chose  between  emission  or  absorbtion  ;  5)  if  at  time  ti  an  OUT  scattering 
is  chosen  the  integration  over  the  final  states  is  performed  analytically,  and  the 
current  k  is  kept  to  continue  the  simulation;  if  an  IN  scattering  is  chosen  an 
initial  state  is  chosen,  at  random  in  the  surface  of  energy  conservation,  and  this 
new  k  is  used  to  continue  the  simulation;  6)  the  previous  two  points  are  repeated 
until  t  =  0  ;  7)  the  initial  distribution  is  evaluated  at  k,(0)  and  multiplied  by 
the  product  of  the  transition  rates  for  the  simulated  path  and  divided  by  the 
probabilities  of  the  choices  selected  in  the  simulation. 

However,  a  change  can  be  introduced  in  the  above  algorithm  giving  a  more 
physical  and  fast  evaluation  of  f.  In  fact  all  terms  containing  any  number 
of  OUT  scattering  events  between  two  given  IN  scattering  can  be  summed 
up  analitically.  Let  us  consider,  for  example,  the  term  containing  two  OUT 
scattering  events  : 

J  *i  J  dk,P0(k(<i),ki)  J  dt2  f  dk2Po(k(t2),k2)/(k(0),0)  (2.46) 

Since  the  integrals  are  symmetric  with  respect  to  the  echange  of  the  time  vari¬ 
ables,  Eq  (2.26)  can  be  reduced  to  : 

/(k(0),°)i  J  dh  j  dk1Pc(k(<l),k1)  j dt,  j  dk2P0(k(t2),k2 
=  /0^[/*1Po(k,f1)]2 

where  P0  (k , 1 1 )  is  the  scattering  probability  integrated  over  all  possible  final 
states. 
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If  we  sum  all  such  contribution  the  result  is  the  exponential  series: 

=  /0exp{-|at1Po(k1,t)}  (2.47) 

Where  now  exp  f  df  iP0(k,  ti)  is  the  probability  of  no  scattering  between  the 
times  0  and  t.  new  term  is  the  ballistic  term  moltiplied  by  the  probability  that 
am  electron  has  not  had  a  scattering  in  the  time  interval  (0 ,t). 

The  sarnie  argument  can  be  applied  to  sum  up  analytically  all  OUT  terms 
between  two  given  IN  events.  The  result  is  that  for  any  pair  of  IN  terms  at  times 
1 1  and  £2  the  estimator  must  be  moltiply  by  the  probability  exp{—  f  <ftP0(k,t)}, 
that  an  electron  is  not  taken  away  from  the  trajectory  by  an  OUT  scattering. 

The  procedure  sketched  above  is  then  to  be  modified  in  that  only  IN  scat¬ 
tering  events  are  to  be  considered  and  the  no-scattering  probabilities  due  to 
out  scattering  added  to  the  estimator.  Such  a  change  results  in  an  enormous 
improvment  of  the  convergence  since  we  move  from  a  series  with  alternate  signs 
to  a  series  of  positive  terms. 

2.2.5. 2  Results 

The  above  method  has  been  applied,  for  testing  purposes,  to  a  simple  model 
semiconductor  based  on  silicon.  A  simple  spherical  parabolic  band  is  used  and 
only  one  optical  phonon  scattering  is  considered.  For  such  a  model  the  integral 
in  equation  (2.47)  can  be  easily  performed  analytically.  The  following  physical 
parameters  have  been  used:  effective  mass=.295,  phonon  temperature=  450° A', 
coupling  constant  0.8  109  V/cm,density=2. 329cm-3  .  Figs  (2.21)  and  (2.22) 
show  the  results. 

Fig  (2.21)  shows  the  distribution  function  obtained  at  different  times  at  the 
initial  application  of  the  electric  field  E=  10KV /cm.  While  for  t  —  0.5ps  the 
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curve  is  complete  with  terms  up  to  fifth  order,  for  t  =  1.5 ps  we  had  to  include, 
for  convergence,  terms  up  to  the  12th  order,  (path  with  up  to  12  scattering 
events  contribute  to  the  distribution  function)  In  Fig  (2.22)  the  distribution 
function  is  shown  as  a  function  of  energy.  A  sampling  of  final  k  values  in  the 
energy  sphere  had  to  be  performed.  Results  of  the  new  method  are  compared 
with  a  traditional  (EMC)  performed  with  2.4  106  particles.  The  curves  in  the 
figure  show  how  many  orders  contribute  to  the  final  shape  of  the  distribution 
function. 

Several  problems  could  be  discussed  about  the  efficiency  of  the  method 
proposed  here.In  particular  it  is  to  be  considered  that  if  the  probabilities  (the 
p,’s  that  appear  in  Eq  (2.22))  are  not  appropriately  chosen  the  advantage  of 
selecting  the  final  k  is  lost  by  the  fact  that  most  of  the  times  the  initial  kj, 
where  the  source  function  /(ki(0),0)  is  to  be  evaluated,  falls  into  regions  where 
f  is  very  small  and  the  variance  would  result  very  large.  A  smart  choice  of 
the  pi’s  for  every  particular  problem  will  greatly  improve  the  efficency  of  the 
method. 

Finally  we  would  like  to  remind  that  the  method  should  be  applicable  also 
to  space  dipendent  problems;  in  particular  it  could  be  very  useful  to  analyse  the 
problem  of  electrons  injected  into  the  semiconductor  material  from  a  metallic 
contact.  In  the  traditional  M.C.  very  long  paths  (and  computations)  correspond 
to  electrons  wandering  around  in  the  metal  before  entering  the  semiconductor 
space.  On  the  contrary,  with  the  present  procedure  only  paths  ending  in  the 
desired  region  are  analysed. 


35 


References 


[2.1]  For  a  recent  review  on  the  subject  see  L.  Reggiani,  Physica  134B,  123 
(1985) 

[2.2]  I.B.  Levinson  and  Ya.  Yasevichyute,  Sov.  Phys.  JETP  35,  991  (1972) 

[2.3]  J.R.  Barker,  Solid  State  Electron.  21,  267  (1978). 

[2.4]  K.K.  Thornber,  Sol.  State  Electron.  21,  259  (1978) 

[2.5]  D.C.  Herbert  and  S.J.  Till,  J.  Phys.  C,  15,  5411  (1982) 

[2.6]  N.  Pottier  and  D.  Calecki  Physica  110A,  471  (1982) 

[2.7]  V.P.  Seminozhenko,  Phys.  Repts.  3,  103  (1982) 

[2.8]  A.C.  Marsh  and  J.C.  Inkson,  J.  Phys.  C:  Solid  State  Phys.  17,  4501  (1984) 

[2.9]  D.J.  Lowe,  J.  Phys.  C:  Solid  State  Phys.  18,  L209  (1985). 

[2.10]  S.K.  Sarker,  J.H.  Davies,  F.S.  Khan  and  J.W.  Wilkins,  Phys.  Rev.  B33, 
7263  (1986) 

[2.11]  F.S.  Khan,  J.H.  Davies  and  J.W.  Wilkins,  Phys.  Rev.  B36,  2578  (1987) 

[2.12]  G.D.  Mahan  in  "Polarons  in  ionic  crystals  and  polar  semiconductors”  Ed. 
J.T.  Devreese,  North-Holland  (Amsterdam,  1972)  p.  554 

[2.13]  G.D.  Mahan  ”Many-particle  Physics”  (Plenum,  New  York,  1981) 

[2.14]  J.R.  Barker,  J.  Phys.  C:  Solid  State  Phys.  6,  2663  (1973) 

[2.15]  F.  Capasso,  T.P.  Pearsall  and  K.K.  Thornber,  IEEE  EDL-2,  295  (1981) 

[2.16]  J.Y.  Tang,  H.  Shichijo,  K.  Hess  and  G.J.  Iafrate,  Journal  de  Physique  Coll. 
-  C7,  42,63  (1981). 

[2.17]  Y.C.  Chang,  D.Z.Y.  Ting,  J.Y.  Tang  and  K.  Hess,  Appl.  Phys.  Lett.  42, 
76  (1983) 

[2.18]  J.Y.  Tang  and  K.  Hess,  J.  Appl.  Phys.  54,  5139  (1983) 

[2.19]  K.  Brennan  and  K.  Hess,  Solid  State  Electron.  27,  347  (1984) 

[2.20]  K.  Brennan  and  K.  Hess,  Phys.  Rev.  B29,  5581  (1984) 


36 


[2.21]  S.D.  Brorson,  D.J.  DiMaria,  M.V.  Fischetti,  F.L.  Pesavento,  P.M.  Solomon 
and  D.W.  Dong,  J.  Appl.  Phys.  58,  1902  (1985) 

[2.22]  M.  Artaki  and  K.  Hess,  Superlattices  and  Microstructures  1,  489  (1985) 

[2.23]  W.  Porod  and  D.K.  Ferry,  Physica  134B,  137  (1985) 

[2.24]  P.  Lugli,  L.  Reggiani  and  C.  Jacoboni,  Superlattices  and  Microstructures 
2,  143  (1986) 

[2.25]  L.  Reggiani,  P.  Lugli  and  A.P.  Jauho,  Phys.  Rev.  B36,  6602  (1987) 

[2.26]  A.P.  Jauho  and  J.W.  Wilkins,  Phys.  Rev.  29B,  1919  (1984). 

[2.27]  D.  Langreth  and  J.  Wilkins,  Phys.  Rev.  6B,  3189  (1972);  D.C.  Langreth, 
in  "Linear  and  Nonlinear  Electron  Transport  in  Solids”,  Eds.  J.T.  Devreese 
and  E.  Van  Doren  (Plenum,  New  York,  1976). 

[2.28]  P.  Lipavsky,  V.  Spicka  and  B.  Velicky,  Phys.  Rev.  34B,  6933  (1986) 

[2.29]  G.D.  Mahan,  Phys.  Rep.  110,  321  (1984);  ibidem  145,  235  (1987) 

[2.30]  A.P.  Jauho  and  L.  Reggiani,  Solid  State  Electron,  to  be  published 

[2.31]  C.  Jacoboni  and  L.  Reggiani,  Rev.  Mod.  Phys.,  55,  645  (1983) 

[2.32]  L.  Reggiani,  "Hot  electron  transport  in  semiconductors" ,  Topics  in  Applied 
Physics  58,  Springer  Verlag  (Heidelberg,  1985) 

[2.33]  J.  Liu  and  L.C.  Chin,  Appl.  Phys.  Lett.  47,  1304  (1985) 

[2.34]  K.  Kim,  B.A.  Mason  and  K.  Hess,  Phys.  Rev.  B36,  6547  (1987) 

[2.35]  T.  Kurosawa,  J.  Phys.  Soc.  Japan,  suppl.21,  424  (1966) 

[2.36]  R.  Brunetti  and  C.  Jacoboni,  Proc.  Int.  Conference  on  Ho  Carriers  in 
semiconductors,  Boston  (1987),  in  press. 

[2.37]  C.  Jacoboni,  P.  Poli,  and  L.  Rota,  Proc.  Int.  Conference  on  Hot  Carriers 
in  semiconductors,  Boston  (1987),  in  press. 

[2.38]  J.M.  Hammersley  and  D.C.  Handscombe,  "Monte  Carlo  Methods”, 
Methuen  and  Co.,  London  (1964) 

[2.39]  J.R.  Barker,  Sol.  State  Electron.  21,  267  (1978) 


37 


Figure  Captions 


Fig.  2.1  -  Joint  spectral  density  accounting  for  ICFE  as  a  function  of  the  kinetic 
j  energy  after  a  scattering  event.  Continuous  and  dashed  lines  refer  respectively 

to  the  exact  and  approximated  expressions  reported  in  the  text  for  an  initial 
:  energy  e,  =  leV.  with  E  =  500 kV/cm,  and  q  antiparallel  to  E. 

i  Fig.  2.2  The  same  as  Fig.  2.1  with  q  parallel  to  E. 

Fig.  2.3  The  same  as  Fig.  2.1  with  E  =  100 kV/cm,  and  q  parallel  to  E. 

!  Fig.  2.4  -  Joint  spectral  density  accounting  for  CB  as  a  function  of  the  kinetic 

energy  after  a  scattering  event.  Continuous  and  dashed  lines  refer  respectively 
to  the  exact  and  approximated  expressions  reported  in  the  text  with  -y2  = 
l.lmeV.and  for  an  initial  energy  £j  =  leV . 

Fig.  2.5  -  The  same  as  Fig.  2.4  with  €,•  =  0.2eK. 

Fig.  2.6  -  The  same  as  Fig.  2.5  with  72  =  O.llmeF. 

Fig.  2.7  -  Distribution  function  of  the  carrier  kinetic  energy  at  E  =  500 kV/cm. 
Dashed  curve  refers  to  a  semiclassicai  (SC)  simulation,  continuous  curve  to 
a  simulation  which  includes  collisional  broadening  (CB)  only,  and  dot-dashed 
curve  to  a  simulation  which  includes  intra-collisional  field  effects  (ICFE)  only. 

Fig.  2.8  -  The  same  as  Fig.  2.7  at  E  =  100 kV/cm. 

Fig.  2.9  -  The  same  as  Fig.  2.7  at  E  =  10 kVjcm 

Fig.  2.10  -  The  same  as  Fig.  2.7  at  E  =  5 kV/cm 
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Fig.  2.11  -  Distribution  function  of  the  difference  between  the  initial  and  final 
kinetic  energy  after  a  scattering  event  (Ac  =  f  j  —  c/)  as  obtained  from  Monte 
Carlo  simulations  for  the  case  E  =  500 kV/cm.  The  vertical  line  at  Ac  = 
hu> o  represents  the  delta  distribution  of  the  semiclassical  case.  The  continuous 
(dashed)  curves  refer  to  CB  and  (ICFE),  respectively. 

Fig.2.12  -  Diagrams  representing  the  second  order  contributions  to  the  density 
matrix.  The  horizontal  axes  represent  the  time  for  the  two  arguments  of  the 
density  matrix  and  arrows  indicate  phonon  absorption  and  emission  processes. 

Fig.  2.13  -  Diagram  representing  the  fourth-order  contribution  to  the  density 
matrix  shown  in  Eq.(2.30). 

Fig.  2.14  -  Absolute  value  of  quantum  corrections  at  the  second  perturbative 
order  compared  with  the  absolute  value  of  the  classical  one-scattering  correction 
for  the  model  semiconductor  (see  text) .  (a)  classical  one-scattering  contribution; 
(b)  quantum  correction  of  the  second  perturbative  order;  (c)the  same  as  in  (b) 
without  ICFE. 

Fig.  2.15  -  Quantum  corrections  at  the  fourth  perturbative  order  compared 
with  the  classical  two-scattering  correction  for  the  model  semiconductor  (see 
text),  (a)  classical  two-scattering  contribution;  (b)  quantum  correction  of  the 
fourth  perturbative  order;  (c)  the  same  as  in  (b)  with  separate  collisions;  (d) 
the  same  as  in  (c)  without  ICFE. 

Fig.  2.16-  Density  matrix  to  zero-order  (a),  to  second-order  (c),  and  to  fourth- 
order  (b)  perturbative  correction  for  the  silicon-like  model  (see  text). 
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Fig..  2.17  -  Classical  distribution  function  of  electrons  as  a  function  of  energy 
at  a  time  t  =  Chips  after  excitation.  The  highest  peak  at  1000/f  is  the  initial 
distribution  at  t  =  0. 

Fig.  2.18  -  Quantum  distribution  function  of  electrons  as  a  function  of  energy 
at  a  time  t  =  O.lps  after  excitation. 

Fig.  2.19  -  Mean  electron  energy  as  a  function  of  time,  obtained  with  quantum 
transport  theory  (continuous  line)  and  classical  theory  (dashed  line). 

Fig.  2.20  -  A  possible  trajectory  of  third  order  with  two  IN  scatterings  and  one 
OUT  scattering. 

Fig.  2.21  -  Electron  distribution  as  a  function  of  wavevector  k  at  different  times 
afetr  the  application  of  the  electric  field  (a:  t=0;  b:  t=0.1  ps;  c:  t=0.5  ps;  d: 
t=l.  ps;  e:  t=1.5  ps). 

Fig.  2.22  -  Electron  distribution  as  a  function  of  energy  obtained  at  time  t=l. 
ps  after  the  application  of  the  electric  field.  Continuous  curves  represent  results 
obtained  by  summing  orders  up  to  10th  (a),  13th  (b),  16th  (c),  19th  (d),  25th 
(e).  E=  105  V/cm. 
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TABLE  2.1 


E 

<(>  " 

Vd 

^maz 

(kV/cm) 

(eV) 

(10  7cm/s) 

(eV) 

0.015 

1.05 

<  1 

SC 

5 

0.018 

0.86 

<  1 

ICFE 

0.015 

1.06 

<  1 

CB 

0.019 

1.06 

<  1 

SC 

10 

0.022 

0.86 

<  1 

ICFE 

0.019 

1.08 

<  1 

CB 

0.13 

1.33 

49 

SC 

100 

0.020 

1.07 

59  ±1 

ICFE 

0.15 

1.38 

110  ±  13% 

CB 

2.5 

1.41 

1200  ±  5% 

SC 

500 

4.7 

1.23 

2200  ±  40% 

ICFE 

' 

3.9 

1.65 

3850  ±  13% 

CB 

Table  2.1  -  Summary  of  the  carrier  mean  kinetic  energy  <  c  >,  drift  velocity  Vd, 
and  the  maximum  kinetic  energy  achievable  by  a  carrier  during  a  simulation 
(■max  corresponding  to  the  various  joint  spectral  density  models  at  different 
electric  field  strengths.  The  uncertainty  of  the  numerical  results  is  within  ±2% 
if  not  stated  otherwise. 


TABLE  2.2 


xf 

<  t  > 

{eV) 

Vi 

(10  7cm/s) 

fmai 

(eV) 

2  Xi 

0.20 

1.07 

59 

4ii 

0.29 

0.98 

121 

8  Xi 

0.54 

0.93 

508 

Table  2.2  -  Dependence  of  the  mean  values  and  of  the  maximum  kinetic  energy 
achievable  by  a  carrier  during  a  simulation  imaz  upon  different  values  of  the 
cut-off  energy  for  the  case  of  ICFE  at  E  =  lOOkV/cm 


3.CORRELATION  FUNCTIONS  OF  HOT-ELECTRONS 

3.1  Correlation  Functions  for  Bulk  Semiconductors 

3.1.1  Introduction 

Since  the  work  of  Green  and  Kubo  ,  the  use  of  the  velocity  autocor¬ 
relation  function  has  been  proven  to  be  a  fundamental  tool  for  the  description 
of  the  kinetic  coefficients  under  linear  response  condition  in  an  applied  external 
.  field.  Near  equilibrium,  the  velocity  autocorrelation  function  represents  the  mi¬ 
croscopic  quantity  which  unifies  the  physical  interpretation  of  both  the  mobility 
H  and  diuusivity  D  (fluctuation-dissipation  theorem).  This  is  no  longer  the  case 
under  far  from  equilibrium  conditions,  when  a  static  electric  field  E  externally 
applied  is  sufficiently  high  to  produce  hot-electron  effects. 

Under  far  from  equilibrium  conditions  it  is  possible  to  introduce  a  set  of 
correlation  functions  which  still  enables  a  unified  interpretation  of  the  transport 
coefficients.  In  the  following  sections  we  will  present  theoretical  developments 
and  results  for  bulk  materials. 

.  * 

3.1.2  Equations  of  Motion  for  a  Closed  Set  of  Correlation  Functions 

Following  the  results  of  Niez  [3.1-3]  we  use  a  first  principle  Heisenberg  pic¬ 
ture  to  write  the  equations  of  motion  for  the  required  quantities,  in  our  case 
energy  and  momentum.  By  a  projection  operator  technique,  the  thermodynam¬ 
ical  motion  is  separated  from  the  fluctuations.  Taking  for  simplicity  a  single 
spherical  and  parabolic  band  model  semiconductor,  it  has  been  proven  [3.3] 
that  four  correlation  functions,  which  couple  velocity  and  energy  fluctuations, 
"are  needed.  These  can  be  written  in  matrix  notation 

( <£v*(0  \  _  (  <  6v6v(t)  >  < 

\^tv(i)  4>',{t) )  \<<5c5v(t)>  < 

The  set  in  Eq.  (3.1)  is  found  to  satisfy  the  following 


as: 


6e6e(t)  >  )  ^ 

closed  system  of  coupled 
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first  order  differential  equations: 


Zitvv  = 
2t4>ve  = 


—ai\4>vv  —  C*i2<t>v, 
—021<t>vv  ~  “22^u* 


Tt^'v  = 
li^tt  - 


-Oil <t>tu  -  «12^cc 
-021^>tu  -  CL22<i>tt 


(3.2) 


where  the  coefficients  a,y  (which  depend  on  the  external  field  but  not  on  time) 
describe  the  microscopic  properties  of  the  physical  system.  These  coefficients 
come  from  the  general  theory  [3.3]  when  the  relaxation  time  approximation  is 
used.  Then,  all  the  slow  components  of  the  system  dynamics,  such  as  the  total 
energy  and  momentum,  are  taken  care  of  by  the  projection  operator  used.  The 
kernels  of  the  Langevin-like  convolution  integrals  obtained  in  this  framework 
are  rapidly  varying  functions  [3.4],  and  the  Oy  are  just  their  time  integrals. 
The  physical  meaning  of  the  diagonal  terms  an  and  022  is  related  to  the  mo¬ 
mentum  and  energy  relaxations  rates.  The  off-diagonal  terms  012  and  a2i  are 
zero  at  equilibrium  (E  =  0)  when,  as  known,  the  relaxation  of  momentum  and 
energy  are  independent.  Under  hot  electron  conditions  these  rates  are  no  longer 
independent  and  the  off-diagonal  terms  describe  the  coupling  between  them.  A 
microscopic  determination  of  the  a coefficients  remains  a  formidable  problem 
which  has  not  been  yet  considered.  The  analytical  solution  of  Eqs.  (3.2)  in 
normalized  form  is  given  by: 


w*  >0 
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4>w(t)  =  i - {&2i[woch(wot)  +  (A  —  c*i  i)sh  (w0t)J  +  [(A  —  an)2  —  Wol5MwoO} 

a2iWo 

^**(0  = - (w0ch(w0t)[(A  -  <*n)  +  a21]s/i(woi)]} 

wo 

e-xt 

4>, „(«)  =  - {a2i[woc/i(wo0  +  (A  -  aii)sft(a'0l)|  +  [(A  -  au}2  -  Wo|sfc(w0«)} 

021^0 

^**(0  = - -{w0ch(w0t)[(A  -  an)  +  a2i]a^i(w0<)]} 

Wo 

(3.3a) 


wg  <  0 


e-xt 

^»v(0  =  t - -{a2i[w0eos(w0t)  i-  (A  -  au)sin(wot)]  +  [(A  -  an)2  -  w2]sm(w0f)} 

a2iwo 

g-Xt 

4vc(t)  =  - (w0co.s(wof)[(A  -  an)  +  a2i]s:n(w0J)]} 

U>0 

e- At 

*«„(<)  =  - - {a2i[w0cos(w0i)  +  (A  -  an)sin(wot)]  +  [(A  -  an)2  +w|]sinfw0t)} 

a2iWo 

e~xt 

4et(t)  =  - -{w0cos(w0f)[(A  -  an)  +  a2i]4m(w0I)i} 

Wo 

(3.36) 

where  A  =  l/2(an  +  a22)  ,  wj  =  |  (an  -  a22)2  +  a12a2i  ,  0  = 

<  SvSe  >/<  6v 2  >;  -7  ■—  <  SvSe  >/<  6e2  >;  a2i  =  <*2i/0  ;  a2i  —  a2i-7.  We 
notice  that  A  and  wo  are  two  frequencies  which  characterize  the  main  features  of 
the  time  evolution  of  the  correlation  functions.  As  a  matter  of  fact  A  ,  which  is 
always  real,  is  responsible  for  a  damping  ,  while  wo,  when  imaginary,  determines 
an  oscillatory  behaviour. 
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3.1.3  Numerical  Results 

To  check  the  reliability  of  the  functional  form  given  by  Eqs.  (3.3),  we  have 
prepared  an  Ensemble  Monte  Carlo  simulation  for  the  two  following  cases  of 
interest. 

(i)  Quasi  elastic  regime:  This  condition  corresponds  to  a  carrier  momentum 
distribution  function  with  the  even  part  dominating  over  the  odd  part.  We  have 
considered  here  scatterings  with  acoustic  and  non-polar  optical  (or  intervalley) 
phonons  at  300  K  for  a  simple  model  semiconductor  (Si-like  with  E  //  <  111  > 
[3.5]). 

(ii)  Streaming  motion  regime:  In  this  condition  the  odd  part  of  the  distri¬ 
bution  strongly  prevails  over  the  even  part.  To  this  purpose,  with  respect  to 
the  previous  case,  the  optical  phonon  coupling  has  been  arbitrarily  increased 
(by  a  factor  5)  at  77  K. 

The  results  for  the  above  regimes  axe  shown  in  Figs.  3.1,2,  where  the  data 
obtained  from  the  simulation  are  found  to  agree  quite  satisfactorily  with  the 
analytical  curves  (we  notice  that  a  similar  agreement  has  been  obtained  at 
various  filed  strengths  from  5  to  200  kV/cm).  In  particular  the  quasi  elastic 
regime  is  characterized  by  A  >  |w0|  ,  while  the  streaming  motion  regime  by 
A  <  |wo|  •  In  the  latter  case,  in  agreement  with  the  expectations  [3.6],  the  re¬ 
lationship  w0  =  (2n/c)E(2mhu>op)  h,uiop  being  the  optical  phonon  energy 

and  m  a  carrier  effective  mass,  has  been  verified  within  the  numerical  uncer¬ 
tainty.  Figure  3.3  shows  the  depedence  with  electric  field  of  the  a,/  coefficients 
for  the  case  of  quasi  elastic  regime. 

Having  succeded  in  deriving  a  set  of  equations  describing  the  evolution  of  the 
correlation  functions  of  the  thermodynamical  variables  (v,e),  we  consider  now, 
in  the  same  theoretical  framework,  the  basic  ideas  which  underlie  the  linear 
response  theory  around  the  steady  state.  We  thus  apply  to  the  system  an 
extra  perturbative  electric  field  E\(t).  The  Hamiltonian  is  now  time  dependent, 
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and  we  focus  on  the  thermodynamical  differences  between  the  steady  state 
values  of  our  variables  and  their  values  induced  by  Ei(t).  By  repeating  for  the 
present  case  the  calculations  that  led  to  Eqs.  (3.3),  it  appears  that  the  equation 
which  drives  the  momentum  of  the  system  gives,  as  it  is  well  known,  the  first 
fluctuation-dissipation  theorem  only  when  E  —  0,  Le.  in  the  near  equilibrium 
situation.  If  E  =  0  one  comes  with  the  so-called  “dressing  term”  Te  [3.5],  which 
characterizes  the  dissipative  nature  of  the  steady  state.  We  have  been  able  to 
evaluate  Te  and  to  give  its  espressions  in  terms  of  the  coefficients  driving  the 
dynamics  of  the  correlation  functions  of  the  thermodynamical  variables  only. 
After  some  delicate  algebra  (which  will  be  published  later),  one  can  write  : 


«4(ftc)r‘),,.]E  («-■>.,< ($;!))}  (M 

where  the  notation  of  Niez  [3.3]  has  been  used.  In  Eq.  (3.4)  the  indices  I 
and  2  correspond  respectively  to  the  velocity  and  the  energy  of  the  electronic 
system;  furthermore  is  the  steady-state  drift  velocity.  The  a  coefficients 
are  the  same  as  in  Eq.  (3.2).  The  dressing  term  (second  term  in  the  squared 
braces)  has  been  evaluated  taking  into  account  that  the  collision  time  rc  is  much 
shorter  than  any  time  describing  the  dynamics  of  the  relevant  variables.  In  Te, 
the  ( P\P)l  correlation  functions  of  energy  and  momentum  at  time  t  =  0  are 
defined  through  the  displaced  Maxwellian  distribution  of  the  problem.  The  7 
coefficients  in  Te  are  directly  related  (through  a  mass  factor  m)  to  the  first 
moments  of  the  collision  kernels  [3.3]  for  their  definition): 


t{*J)  ~ 


Hdt  t[G0(t)(l-Uz)  PitL  ;  (1  -  n,)  P)t L]t 
Jo 


(3.5) 


Thus  from  Eq.  (3.4)  it  is  possible  to  give  the  complex  differential  mobility  n'  (w) 
of  the  system  through  information  coming  from  the  correlation  functions  of  the 


thermodynamical  variables  only.  While  an  explicit  formula  for  p'  from  Eq. 
(3.4)  is  currently  in  progress,  its  expression  when  the  dressing  term  is  neglected 
(  i.e.  under  warm  electron  conditions)  is  given  by: 


„  ,  e  r  2A  -  an  -  iu  ^ 


3.1.4  Microscopic  Expression  of  the  Noise  Temperature 

Under  far  from  equilibrium  conditions  due  to  the  presence  of  a  uniform 
electric  field  E,  the  noise  associated  with  velocity  fluctuations  is  commonly 
described  by  the  noise  equivalent  temperature  Tn,  a  physical  quantity  which 
can  be  easily  measured  [3.G],  When  two  particle  interactions  are  neglected  and 
the  system  is  electrically  stable  (i.e.  n'  >  0  ,  n'  being  the  differential  mobility), 
a  generalized  Einstein  relationship  was  proven  to  hold  for  Tn  [3.7]: 


(3.7) 


Here  e  is  the  unit  charge,  Kb  the  Boltzmann  constant  and  D  the  diffusion 
coefficient. 

The  objective  of  this  section  is  to  provide  for  the  first  time  a  microscopic 
expression  for  Tn.  The  essence  of  our  result  relies  on  the  proof  that  both  D  and 
n',  and  therefore  T„,  can  be  defined  from  the  same  set  of  correlation  functions 
^escribing  the  system  under  far  from  equilibrium  conditions. 

To  this  end  we  shall  consider  Eq.  (3.6)  in  the  limit  w  — »  0  and  obtain: 


2A  -  an . 
A*  -  w*  1 


(3.8) 


We  recall  that  Eq.  (3.8)  is  still  an  approximation  since  it  neglects  the  so-called 
“dressing  term”  of  Ref.  [3.3],  which  comes  from  the  dissipative  nature  of  the 
steady-state.  This  limits  the  reliability  of  the  present  theory  to  fields  below  20 
kV/cm. 
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From  the  definition  of  D  in  terms  of  4>vv  [3.6]  we  obtain  the  desired  defini¬ 
tion  of  the  noise  temperature  in  terms  of  the  ax:  coefficients  as: 

r*-5?f  (30) 

Figure  3.4  shows  the  comparison  between  this  theory  and  experiments  [3.8]  of 
the  excess  noise  temperature  for  the  case  of  electrons  in  Si.  The  agreement 
found  is  good  and  supports  the  present  interpretation  of  the  noise  temperature 
in  terms  of  correlation  functions. 

In  summary,  explicit  formulae  have  been  obtained  for  the  time  evolution  of 
a  closed  system  of  correlation  functions  which  generalized  the  linear  response 
theory  to  the  far  from  equilibrium  case  (hot  electrons).  This  theory  provides  a 
unified  microscopic  interpretation  of  both  diffusion  (given  by  the  Fourier  trans¬ 
form  of  <j>uv{t))  and  mobility.  At  zero  field  our  formulation  coincides  with  the 
first  form  of  the  fluctuation-dissipation  theorem.  Therefore,  the  present  theory 
can  be  viewed  as  a  generalization  of  the  Kubo  formalism  under  far  from  equi¬ 
librium  conditions.  As  such,  we  have  proven  that  this  formalism  consistently 
interprets  numerical  simulations,  as  obtained  through  standard  Monte  Carlo 
procedure,  for  a  variety  of  physical  situations  and  experimental  results. 

3.2  Correlation  Functions  for  Quantum  Wells 

3.2.1  Introduction 

A  theoretical  analysis  of  velocity  fluctuations  can  yield  relevant  information 
on  the  microscopic  interpretation  of  transport  coefficients  as  well  as  on  the 
detailed  features  of  the  scattering  sources.  This  analysis  becomes  even  more 
important  if  submicrometer  devices  are  considered,  where  ultrafast  transport 
processes  are  usually  involved  and  a  deeper  insight  into  the  physics  of  transport 
phenomena  is  required.  Furthermore  fluctuations  can  play  an  important  role  in 
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the  design  and  characterization  of  the  device  itself. 

Several  papers  have  appeared  in  the  literature  on  this  subject  for  bulk 
materials,  while  very  scarce  theoretical  investigations  have  been  done  so  far 
in  this  field  for  2D  systems,  and  few  experimental  data  is  available  for  noise 
properties  in  quantum  wells. 

The  theoretical  problem  can  be  approached  through  a  unified  theory  for¬ 
mulated  in  terms  of  the  autocorrelation  function,  a  quantity  directly  related  to 
diffusivity  and  noise.  This  method  cam  be  used  to  describe  both  steady-state 
(Sec.3.2.2)  and  transient  (Sec. 3.2.3)  situations,  and  also  to  analyse  the  different 
contributions  due  to  the  different  physical  sources  of  fluctuations  which  arise  in 
the  presence  of  an  applied  electric  field. 

In  order  to  obtain  results  for  realistic  structures,  a  Monte  Carlo  simulation 
of  the  quantum  well  has  been  used,  that  allows  to  obtain  an  exact  solution  of  the 
transport  problem  in  hot  electron  conditions  (i.e.  out  of  equilibrium,  in  presence 
of  high  electric  fields),  where  analitical  techniques  cannot  be  succesfully  applied 
without  introducing  severe  approximations.  Furthermore,  a  direct  simulation 
of  the  dynamics  of  charge  carriers  inside  the  crystal  enables  us  to  extract  any 
required  physical  information  while  the  solution  of  the  transport  equation  is 
being  built  up  and  the  simulation  can  be  easily  modeled  in  order  to  reproduce 
particular  experimental  conditions. 

The  physical  system  and  the  details  of  the  Monte  Carlo  procedure  used 
for  the  simulation  of  the  2D  electron  gas  in  the  quantum  well  are  presented  in 
Sec.3.2.4.  Results  for  autocorrelation  function,  diffusivity,  and  noise  at  room 
temperature  are  discussed  in  Sec.  3.2.5  and  3.2.6  for  stationary  conditions  and 
transient  conditions,  respectively,  together  with  a  comparison  between  2D  and 
3D  data. 

Few  syntetic  conclusive  statements  are  given  in  Sec. 3. 2. 7  together  with 
possible  developments  of  this  research  for  the  future. 
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3.2.2  Autocorrelation  Function,  Diffusion  and  Noise  in  Stationary  Con¬ 
ditions 

Let  us  consider  an  ensemble  of  electrons  in  a  crystal  semiconductor  subject 
to  an  external  electric  field  E  and  to  the  action  of  scattering  agents  (phonons, im¬ 
purities,  etc).  If  the  electrons  have  any  nonuniform  distribution  in  space  the 
phenomenon  of  diffusion  occurs,  tending  to  make  the  concentration  uniform 
through  the  spreading  out  of  the  carriers. 

Diffusion  is  described,  at  a  phenomenological  level,  by  the  equation  [3.9]: 

Ji  =  e{n(r)»rf,  -  (3.10) 

In  the  above  equation  e  is  the  electron  charge,  r  is  the  space  position  with 
components  z,,  n(r)  and  J  the  particle  density  and  current  density,  respectively, 
Dij  is  the  diffusion  coefficient  tensor  (ij=l,2,3,  the  sum  over  repeated  indices 
is  implied),  Vj  the  drift  velocity  of  the  carriers  in  absence  of  diffusion. 

If  E  is  applied  along  a  high  symmetry  direction  of  a  cubic  crystal,  then 
Dij  reduces  to  a  diagonal  form,  with  a  longitudinal  component  Di  and  two 
transverse  components  Dt.  In  the  following  we  will  analize  diffusivity  and  fluc¬ 
tuations  along  the  field  direction;  consequently  we  will  use  a  simplified  scalar 
notation  where  all  the  vector  quantities  are  substituted  with  their  longitudinal 
components. 

For  vanishing  small  electric-field  strengths,  diffusivity  D  and  mobility  n 
are  field  independent  and  verify  the  Einstein  relation: 


D  hKbT 
e 


(3.11) 


where  T  is  the  crystal  temperature  and  Kb  the  Boltzmann  constant. 

At  high  fields  the  Einstein  relation  fails  and  the  study  of  diffusion  is  gen¬ 
erally  performed  through  the  introduction  of  a  field-dependent  D  [3.10]. 
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If  the  concentration  gradients  are  small,  in  absence  of  carrier-carrier  inter¬ 
actions  D(E)  can  be  obtained  from  the  following  equation: 

=2D  (3.12) 

where  x  is  the  displacement  along  the  field  direction  and  brackets  represent 
ensemble  average.  The  quantity  on  the  right  hand  side  is  the  second  central 
moment  of  the  distribution  n(r)  (SCM);  eq.(3.12)  is  valid  at  times  longer  than 
both  the  transient  transport  time  and  the  time  necessary  for  setting  up  the 
correct  space-velocity  correlations  at  the  basis  of  diffusivity  effects. 

The  diffusion  phenomenon  is  strictly  related  to  velocity  fluctuations  and 
noise.  The  mathematical  quantity  that  describes  the  common  origin  of  diffusion 
and  noise  is  in  fact  the  autocorrelation  function  of  velocity  fluctuations,  which 
carries  the  information  on  how  large  the  fluctuations  are  and  how  they  decay 
in  time: 

C(t)  =  (6v(t)6v(t  +  t))  (3.13) 

(the  mean  value  in  steady  state  conditions  is  independent  of  r).  C[t)  is  related 
to  the  diffusion  D  through  [3.9] : 


f°° 

D=  dtC(t ) 

Jo 


Thus  D  can  be  evaluated  from  C(t),  which  is  of  interest  in  itself  since  it  gives 
important  physical  information  on  the  time  evolution  of  the  dynamics  of  the 


Finally,  we  introduce  the  noise  spectrum  S„(w): 


rT  5 

5„(w)  =f«mr_00(  1/  6v(t)e'utdt\ 

Jo 


Another  well-known  relation  between  C(t)  and  the  noise  spectrum  is  given  by 
the  Wiener-Kintchine  theorem  [3.9): 


Sv(w)  ~-2  f  dtC(t)e 
Jo 
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From  this  last  equation  and  eq.(3.14)  extended  to  non-zero  frequencies,  we  have: 

D{u>)  =  (3.17) 

Different  physical  sources  contribute  to  fluctuations  and  diffusion  of  carriers  in 
quantum  wells.  As  it  happens  in  bulk  materials,  fluctuations  in  carrier  momen¬ 
tum  produce  the  so  called  r  thermal”  velocity  fluctuations  [3.9] ;  energy  fluctua¬ 
tions  are  associated  with  "convective”  velocity  fluctuations  and  noise  [3.11],  and 
valley  fluctuations  bring  about  "intervalley”  velocity  fluctuations  [3.12],  when, 
as  in  our  case,  two  or  more  non  equivalent  valleys  exist.  In  the  case  of  quantum 
wells,  however,  the  existence  of  several  subbands  with  different  average  veloci¬ 
ties  implies  the  appearance  of  a  new  source  of  velocity  fluctuations  and  noise, 
that  is  produced  by  the  fluctuation  of  the  subband  occupied  by  the  electrons 
during  their  motion  in  the  quantum  well.  In  the  following  we  will  refer  to  these 
fluctuations  as  to  "intersubband”  velocity  fluctuations. 

Following  a  decomposition  procedure  aiready  applied  to  bulk  structures 
[3.13,3.14[  let  us  consider  an  electron  that,  at  time  t,  is  in  a  subband  of  type  B(t) 
and  in  a  valley  of  type  V(t)  (the  indeces  V  and  B  depend  on  time  because  carriers 
during  their  motion  change  both  valley  and  subband  because  of  scattering);  let 
also  u,vb(0  be  the  mean  velocity  of  electrons  in  valley  V  and  subband  B  with 
energy  between  (  and  f  +  de. 

The  instantaneous  velocity  of  each  electron  v(t)  can  then  be  written  as  the 
drift  velocity  plus  a  number  of  fluctuating  terms: 

v{t)  —  Vd  +  [ufl(0  —  V<i]  +  (vyfl(f)  -  Vfi(t)]  +  [v£vs(0  —  VvcMl  +  [v(0  -  w«VbW] 

=  Vi  +  Svff(t)  +  6w[t)  +  6vt(t)  +  6vK(t),  (3.18) 

Svg(t)  and  6vy(t)  are  the  fluctuations  associated  with  the  drift  velocity  of  the 
subband  and  the  valley  in  which  the  electron  is  at  time  t, respectively;  6vc(t) 
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is  the  velocity  fluctuation  associated  to  the  fluctuation  of  electron  energy,  and 
6vK(t)  is  the  velocity  fluctuation  associated  with  the  fluctuation  of  the  electron 
momentum. 

By  using  the  definition  in  Eq.(3.13),  the  steady-state  autocorrelation  func-1 
tion  becomes 

C(t)  =  +  0>  =  E  (3-19) 

ij  ij 

where 

<?«(*)  =  <MO*V(*'+0>.  (3.20) 

and  i,j  =  B,  V,  e,  k.  The  total  autocorrelation  function  of  velocity  fluctuations 
contains  four  “  diagonal”  contributions  C,,(t)  in  Eq.(3.19),  at  the  origin  of 
interband  ( Cbb ),  intervalley  (CVv),  convective  (Cte),  and  thermal  (C**)  noise, 
respectively.  In  general,  however,  off-diagonal  terms  C{j  also  contribute  to  the 
autocorrelation  of  velocity  fluctuations  [3.13,3.14],  Only  when  the  relaxation 
times  of  the  various  fluctuating  terms  have  well  different  values,  in  calculating 
the  “off-diagonal”  terms  one  of  the  two  fluctuations  can  be  assumed  as  constant, 
while  the  other  fluctuation  averages  to  zero. 

Owing  to  the  linearity  of  Eqs.(3.14)  and  (3.16)  we  can  also  associate  specific 
terms  contributing  to  the  autocorrelation  function  with  corresponding  terms 
contributing  to  dinusivity  and  noise,  thus  making  explicit  their  physical  origins. 

3.2.3  Autocorrelation  Function  and  Diffusion  in  Transient  Conditions 

The  diffusion  process  of  a  carrier  ensemble  comes  from  the  particle  space- 
velocity  correlations  which  arise  during  the  evolution  in  time  of  the  system. 
Starting  from  an  initial  condition  in  which  the  particle  positions  and  velocities 
are  totally  uncorrelated,  the  process  which  occurs  during  the  time  necessary  for 
setting  up  the  correlations  will  be  defined  as  correlation  transient.  Furthermore, 
when  a  high  electric  field  is  applied  at  a  certain  time  to  the  electron  ensemble, 


the  transport  process  itself  must  pass  through  a  transient  region  which  is  nec¬ 
essary  for  attaining  the  stationary  distribution  /( k)  in  k  space.  This  process  is 
the  transport  transient. 

The  definition  of  the  transient  diffusion  coefficient  has  been  given  by  a 
generalization  of  Eq.  (3.12)  to  arbitrary  small  times  [3.15-3.17]: 

=  (3-21) 

where  z(t)  is  the  space  position  of  a  carrier  at  time  t  alon£  the  z  direction 
parallel  to  Vj. 

This  generalization  can  be  put  in  an  equivalent  form  in  terms  of  the  two- 
time  transient  autocorrelation  function.  If  there  is  no  correlation  between  the 
initial  positions  and  velocities  of  the  particles  we  have  [3.16]: 

Z>(t)  =  /  drCt(r),  (3.22) 

Jo 

with 

Ct(r )  =  {6v(t)6v(t  -  r)),  0  <  r  <t  (3.23) 

Eq.(3.22)  reduces  to  Eq.(3.14)  in  steady-state  conditions  (t  — *  co). 

By  comparing  these  two  equations, we  see  that  in  transient  cases  (i)  the  inte¬ 
gration  interval  ranges  from  t=0  (initial  conditions)  up  to  t  (observation  time); 
(ii)  the  autocorrelation  function  to  be  integrated  in  the  transient  analysis  is  not 
time  independent;  it  is  given  by  the  specific  ensemble  average  at  a  particular 
time,  and  its  shape  provides  information  about  the  transport  transient. 

3.2.4.  The  Physical  System  and  The  Monte  Carlo  Procedure 
Electrons  are  simulated  in  a  square  well  representing  the  effective  ID  po¬ 
tential  arising  from  the  band  offset  between  GaAs  and  AizGaAs  (x=.23).  The 
solutions  of  the  wave  equation  for  this  potential  give  rise  to  a  series  of  2D  sub¬ 
bands  which  are  used  in  calculating  the  scattering  rates  for  electronic  motion 
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parallel  to  the  well,  modeled  using  an  Ensemble  Monte  Carlo  simulation.  We 
treat  both  intra-  and  intersubband  scattering  of  the  2D  electrons  by  bulk  lon¬ 
gitudinal  optical  phonons;  intervalley  scattering  to  the  satellite  L- valleys  is  also 
included.  Details  of  the  physical  model  and  of  the  Monte  Carlo  code  for  the 
quantum  well  can  be  found  in  Refs. 3. 18  and  3.19. 

From  the  time  evolution  of  the  electron  ensemble,  we  calculate  the  sec¬ 
ond  central  moment  of  the  carrier  displacement  as  a  function  of  time  said  the 
transient  autocorrelation  function  Ct(r)  [3.14],  The  transient  diffusion  coeffi¬ 
cient  can  be  evaluated  from  them  in  two  independent  ways  using  eq.(3.21)  and 
eq.(3.22). 

The  analysis  of  the  various  contributions  to  the  stationary  autocorrelation 
function  requires  a  previous  MC  evaluation  of  the  mean  velocities  in  eq.(3.18). 
The  stationary  diffusion  coefficient  is  determined  both  from  the  SCM  following 
eq.(3.12)  and  from  the  ACF  using  eq.(3.14).  The  stationary  noise  spectral 
density  is  evaluated  through  eq.(3.16). 

In  order  to  compare  2D  results  with  3D  results  an  Ensemble  Monte  Carlo 
program  for  bulk  GaAs  has  also  been  used  [3.14].  The  physical  model  for 
GaAs  includes  the  same  intravalley  and  intervalley  scattering  sources  as  the 
2D  program,  and  the  input  parameters  of  the  material  have  been  consistently 
chosen. 

3.2.5  Results  for  Stationary  Conditions 

Results  have  been  obtained  at  300  K  for  a  100  A  well  for  different  applied 
electric  fields.  Fig.3.4  and  Fig. 3.5  report  the  longitudinal  diffusion  coefficient 
and  the  drift  velocity  of  electrons  as  functions  of  field  strength  for  both  the  2D 
system  and  bulk  GaAs.  The  absence  of  dissipative  scattering  mechanisms  below 
the  optical  phonon  energy  does  not  permit  the  simulation  of  ohmic  conditions. 
However  the  extrapolation  to  the  low-field  limit  of  the  data  for  vj  and  D  satisfies 
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the  Einstein  relation  within  the  Monte  Carlo  accuracy. 

The  ohmic  mobility  in  the  two  cases  is  not  signiBcantly  different  because 
the  2D  polar  optical  scattering  rate  is  very  close  to  the  same  curve  for  3D. 
Consequently  the  equilibrium  D  is  also  close  for  2D  and  3D. 

Negative  differential  mobility  is  present  in  both  curves,  and  the  threshold 
field  for  electron  transfer  to  the  upper  valley  and  bands  in  2D  is  loweT  than 
in  3D.  At  room  temperature,  for  fields  lower  than  the  threshold,  the  dominant 
scattering  is  essentially  given  by  intravalley  polar  optical  interaction.  When 
the  electron  energy  is  high  enough  to  allow  transfer  to  upper  valleys  and  upper 
subbands  mobility  decreases  with  field. 

The  diffusion  coefficient  is  found  to  decrease  monotonically  in  both  cases 
at  increasing  fields.  At  low  and  intermediate  fields  D  is  larger  in  2D  than  in  3D 
because  carrier  random  velocities  are  larger  in  the  quantum  well  than  in  3D.  At 
higher  fields  the  difference  between  2D  and  3D  is  reduced  and  finally  it  disap¬ 
pears  when  the  electron  energy  is  high  enough  to  guarantee  full  randomization 
of  electron  paths  in  k  space. 

Our  results  differ  from  the  theoretical  data  of  van  Rheenen  et  al.  [3.20,3.21], 
who  found  large  difference  between  2D  and  3D  results  for  both  drift  velocity  and 
diffusivity  vs.  field,  while  the  present  Monte  Carlo  data  seem  to  show  better 
agreement  with  experiments  [3.21],  even  though  a  direct  comparison  would 
require  the  knowledge  of  accurate  values  of  both  experimental  and  theoretical 
low-field  diffusivity. 

As  a  further  confirmation  of  the  above  interpretation,  Fig.  3.7  shows  the 
normalized  ACF  of  velocity  fluctuations  as  a  function  of  time  for  E=2,6,10 
kV/cm  for  the  2D  case.  A  negative  tail  is  present  at  6  and  10  kV/cm  which  is 
due  to  the  dynamical  effect  of  the  transfer  back  and  forth  from  upper  valleys. 
Furthermore  the  time  necessary  to  cancel  the  correlations  between  velocity  fluc¬ 
tuations  is  larger  at  fields  close  to  the  threshold  field,  and  then  it  decreases  at 
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increasing  fields. 

A  comparison  of  the  data  at  10  kV/cm  with  the  3D  data  at  the  same  field 
reported  in  Ref.3.14  shows  that  the  decay  time  of  the  ACF  is  slightly  shorter 
in  the  2D  case.  The  presence  of  different  subbands  in  fact  tends  to  destroy 
the  strong  correlation  which  sets  in  at  fields  above  threshold  for  upper-valley 
transfer  between  large  positive  k  values  (before  intervalley  transfer),  and  large 
negative  k  values  (after  intervalley  transfer)  [3.22]. 

Fig.  3.8  shows  the  ACF  of  velocity  fluctuations  and  the  different  diagonal 
contributions,  as  analysed  in  Sec.  2.2.2,  for  E=7  kV/cm.  Within  this  analysis 
it  is  seen  that  the  thermal  fluctuations  are  in  this  conditions,  much  larger  than 
the  other  fluctuations,  like  it  happens  in  bulk  GaAs  [3.14],  owing  to  the  high 
electron  energy,  and  they  practically  determine  the  value  of  the  autocorrelation 
function  at  short  times  . 

The  diagonal  terms  Cet  and  Cyv  are  smaller  for  an  order  of  magnitude. 
They  are  almost  equal  in  value  at  short  times,  but  the  intervalley  term  slowly 
decays  to  zero  in  about  3  ps,  while  the  convective  term  vanishes  in  less  than 
1  ps.  The  intervalley  term  is  responsible  for  the  long-time  tail  of  the  total 
autocorrelation  function. 

We  can  notice  that  the  decay  time  of  Ckk  is  longer  than  the  decay  time  of 
Cet;  this  fact  is  related  again  to  the  strong  correlation  between  k  values  above 
threshold  discussed  above,  that  is  peculiar  of  polar  materials. 

The  new  term  Cgg  not  present  in  3D  system  is  negligible  in  these  con¬ 
ditions  because  the  drift  velocities  in  different  subbands  are  very  close  to  each 
other.  We  expect  to  see  an  appreciable  contribution  only  when  the  ID  levels 
are  well  separated  in  energy.  Finally  we  can  say  that  the  overall  contribution 
of  off-diagonal  terms  Cij  is  also  negligible  due  to  the  difference  in  characteristic 
decay  times  associated  to  the  different  sources  of  fluctuation  (0.6  ps  for  Cet,  1.8 
ps  for  Ckk,  2.5  ps  for  Cyy). 
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Fig.3.9  shows  the  power  spectra!  density  of  velocity  fluctuations  as  a  func¬ 
tion  of  frequency  at  E=2,5,  and  10  k V /cm.  This  quantity  exibits  the  same 
qualitative  behaviour  of  the  3D  case  [3.14].  The  bump  at  frequencies  around 
103  GHz  is  due  to  the  presence  of  the  negative  part  in  the  ACF.  The  maximum 
is  shifted  towards  higher  frequencies  as  the  field  increases  because  the  negative 
correlations  appear  at  shorter  times  the  higher  the  field. 

3.2.6  Results  for  Transient  Conditions 

Fig.3.10  shows  the  transient  ACF  Ct(r)  as  evaluated  from  the  Monte  Carlo 
simulation  using  eq.(3.23)  for  E=7  kV/cm  at  different  times  t.  The  evolution 
in  time  of  the  shape  of  this  curve  allows  us  to  verify  directly  how  the  scattering 
agents  act  during  the  very  transient  evolution  of  the  system. 

The  first  curve  at  t=  0.5  ps  shows  a  very  large  negative  part  at  r  between 
.3  and  .5  ps.  At  these  short  times  the  strong  negative  correlation  of  the  values 
of  electron  wavevector  due  to  the  cycling  motion  in  k  space  between  central 
and  upper  valley  is  setting  in.  At  longer  times  it  is  smoothed  down  because 
electrons  progressively  visit  larger  zones  of  k  space. 

At  intermediate  observation  times  (t=2,3  ps)  this  correlation  becomes,  at 
the  longest  r’s,  even  positive  before  staoilizing  approximately  over  the  steady- 
state  value  at  t=4  ps. 

The  initial  values  of  the  transient  ACF  at  very  low  r’s  first  increase  because 
of  the  initial  overshoot  in  carrier  velocity  at  observation  times  of  .5  and  1  ps, 
and  then  they  decrease  monotonically  towards  the  steady-  state  value. 

Fig.3.11  shows  the  transient  diffusion  coefficient  as  a  function  of  time  for 
the  same  fields  of  Fig.  3.7.  For  the  lowest  field  (E- -2  kV/cm)  below  threshold 
for  intervalley  transfer  D  reaches  monotonically  the  steady-state  value.  At  5 
and  10  kV/cm  D  first  increases  rapidly,  then  it  decreases  to  values  lower  than 
the  steady-state  value,  and  finally  it  reaches  the  stationary  value.  In  SD  GaAs 
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D  is  found  to  be  ever,  negative  because  of  the  same  effect,  i.e.  the  electron  cloud 
after  a.  first  initial  spread,  for  z.  short  time  shrinks,  and  then  it  is  blown  around 
again.  This  peculiar  behaviour  of  polar  materials  is  related  to  the  decrease  in 
velocity  during  the  transient,  when  the  electron  cloud  begins  the  transfer  to 
upper  valleys  and  bands. 

In  2D  systems  the  presence  of  upper  3ubbands  seems  to  smooth  the  strong 
negative  correlation  of  the  dynamical  picture  described  above,  and  D  is  always 
found  positive. 

Finally  Fig.3.12  shows  the  transient  evolution  of  the  diffusion  coefficient  for 
two  different  carrier  initial  conditions.  Cur  ve  (a)  describes  carriers  starting  from 
an  equilibrium  thermal  distribution  in  the  central  valley  of  the  lower  subband, 
while  curve  (b)  corresponds  to  an  initial  monoenergetic  carrier  distribution  in 
the  upper  valley  of  the  second  subband,  with  energy  much  larger  than  both 
equilibrium  and  steady-state  mean  energy,  that  may  roughly  describe  the  initial 
condition  of  optically  excited  carriers. 

From  the  comparison  of  the  two  cases  it  can  be  seen  that  the  overshoot 
of  D  is  strongly  enhanced  in  the  second  case  (for  about  a  factor  3),  but  the 
undershoot  is  not  strongly  changed  in  value,  even  though  the  minimum  is  shifted 
at  longer  times.  Even  in  this  case  transient  D  is  never  found  negative. 

In  summary,  a  theoretical  analysis  of  velocity  fluctuations  in  GaAs-AlGaAs 
quantum  wells  both  in  steady-state  and  in  transient  conditions  has  been  per¬ 
formed.  Different  contributions  to  velocity  fluctuations  and  diffusion  are  anal¬ 
ysed  separately  in  a  2D  system  and,  in  particular,  the  contribution  of  inter¬ 
subband  fluctuations,  a  new  source  not  present  in  bulk  system,  is  estimated. 
This  study  allows  to  acquire  a  deeper  insight  into  the  microscopic  details  of  the 
transport  picture  in  these  structures. 

A  theoretical  analysis  of  the  transient  properties  of  fluctuations  and  dif- 
fusivity  of  electrons  in  quantum  wells  is  presented  for  the  first  time  .  Results 


show  that  the  transient  features  are  strongly  dependent  on  initial  conditions, 
and  that,  for  initial  hot  carrier  distributions,  D  can  be  three  times  larger  than 
for  the  case  of  equilibrium  distributions.  These  results  are  important  particu¬ 
larly  in  connection  with  the  realization  of  submicrometer  devices,  for  which  the 
transit  time  can  be  comparable  with  the  time  required  to  attain  steady-state 
conditions,  and  the  distance  travelled  by  carriers  during  the  transient  cam  be  a 
significative  portion  of  the  device  length. 

This  amalysis  will  be  extended  in  the  future  to  quantum  wells  of  different 
depths,  and  to  more  realistic  models  including  real-space  tramsfer,  electron- 
electron  and  electron-hole  interactions. 
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Figure  Captions 


Fig.  3.1  -  Normalized  correlation  functions  as  a  function  of  time  for  the  case 
of  electrons  in  Si  at  T  =  300  K  and  E  =  10  kV /cm.  Continuous  and  dashed 
curves  refer  to  present  theory  and  Monte  Carlo  calculations,  respectively. 

Fig.  3.2  -  The  same  as  Fig.  3.1  for  the  case  of  “streaming  motion  regime” 

Fig.  3.3  -  Set  of  cn,  coefficients  as  a  function  of  the  electric  Euld  for  the  case  of 
quasi  elastic  regime. 

Fig.  3.4  -  Excess  no'.se  temperature  as  a  function  of  the  electric  field  in  Si  at  T 
=  300  K.  The  points  refer  to  experiments  of  Ref.  [3.8j  and  the  curve  to  present 
theory.  The  error  bars  indicate  the  uncertainty  of  the  calculations. 

Fig.  3.5  -  Longitudinal  diffusion  coefficient  as  a  function  of  field  strength  for 
2D  and  3D  cases.  Horizontal  lines  indicate  equilibrium  values  for  the  two  cases. 

Fig.  3.6  -  Drift  velocity  as  a  function  of  field  strength  for  bulk  GaAs  and  the 
2D  quantum  well  at  room  temperature. 

Fig.  3.7  -  2D  autocorrelation  function  of  velocity  fluctuations  as  a  function  of 
time  for  E=2,6,  and  10  kV/cm. 

Fig.  3.8  -  Autocorrelation  function  of  velocity  fluctuations  at  7  kV/cm  and  its 
diagonal  associated  terms  as  functions  of  time  for  the  2D  case  (see  text)  . 

Fig.  3.9  -  Power  spectral  density  of  velocity  fluctuations  as  a  function  of  fre¬ 
quency  at  room  temperature  for  E=2,5  and  10  kV/cm  for  the  2D  case. 

Fig.  3.10  -  Transient  autocorrelation  function  of  velocity  fluctuations  for  the 
2D  system  as  a  function  of  r  evaluated  at  different  times  t  indicated  by  the 
numbers  on  each  curve. 

Fig.  3.11  -  Transient  diffusion  coefficient  as  a  function  of  time  for  the  2D  system 
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at  the  same  fields  shown  in  Fig. 3.  Electrons  are  at  t=0  in  equilibrium  conditions. 

Fig.  3.12  -  Transient  diffusion  coefficient  as  a  function  of  time  for  the  2D  system, 
a)  refers  to  equilibrium  initial  conditions,  i>)  refers  to  or.  ir.it i 1,1  monoenergetic 
carrier  distribution  at  103  K. 
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4.  MONTE  CARLO  STUDIES  OF  NON-EQUILIBRIUM  PHONON 
EFFECTS 


4.1  Introduction 

In  the  presence  of  strong  external  perturbations  (laser  excitation,  applied 
electric  fields),  the  energy  transfered  to  a  semiconductors  drives  the  carrier 
system  out  of  equilibrium.  Such  a  situation  is  usually  referred  to  as  hot  car¬ 
rier  condition.  If  the  main  dissipation  channel  for  the  carriers  is  via  phonon 
emission,  then  a  non  equilibrium  (hot)  phonon  population  can  be  found  as  a 
result  of  the  energy  transfer  to  the  lattice.  The  presence  of  phonon  amplifi¬ 
cations  will  ultimately  depend  on  the  rate  at  which  carriers  supply  energy  to 
the  phonons  compared  to  the  rate  at  which  the  phonons  dissipate  their  excess 
energy  to  the  thermal  bath.  The  flux  of  energy  in  and  out  of  the  carrier  system 
is  schematically  shown  in  Fig.  4.1.  In  general,  both  type  of  carriers  (electrons 
and  holes)  can  be  present,  and  channels  for  energy  exchange  are  provided  by 
their  mutual  interaction,  as  well  as  by  their  interaction  with  the  lattice.  While 
Raman  measurements  have  been  able  to  detect  hot  phonon  distributions  for  the 
case  of  picosecond  and  subpicosecond  excitations,  only  indirect  and  non  con¬ 
clusive  evidences  exist  today  on  the  possibility  for  the  phonon  disturbaces  to 
feed  back  into  the  carrier  systems  and  modify  the  overall  energy  transfer  pro¬ 
cess.  As  the  framework  of  semiclassical  transport  theory  is  well  suited  for  the 
description  of  the  transient  dynamics  of  coupled  nonequilibrium  carrier-phonon 
systems,  almost  all  theoretical  work  on  hot  phonons  has  hitherto  considered  the 
Boltzmann  equation  for  the  phonons  with  simplifying  assumptions  about  the 
functional  shape  of  the  carrier  distributions.  The  purpose  of  the  present  contri¬ 
bution  is  to  circumvent  such  assumptions  by  directly  solving  the  coupled  Boltz¬ 
mann  equations  for  phonons  and  carriers  by  a  Monte  Carlo  simulation.  The 
advantages  of  the  Monte  Carlo  technique  are  that  it  provides  a  very  accurate 
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microscopic  description  of  the  physical  processes,  and  it  does  not  require  any 
assumption  on  the  electron  nor  on  the  phonon  distribution  functions.  Further¬ 
more,  a  direct  evaluation  of  the  characteristic  times  for  the  various  scattering 
mechanisms  involved  in  the  interaction  process  is  possible.  The  dynamics  of  the 
carrier-phonon  system  can  be  extremely  complex.  For  instance,  in  the  case  of 
laser  excitation  in  GaAs,  six  types  of  interacting  carriers  ( T,L  and  X  electrons, 
light,  heavy,  and  split-off  holes)  might  be  present.  Depending  on  the  range 
of  excited  densities  and  on  the  lattice  temperature,  different  scattering  mecha¬ 
nisms  can  be  of  importance.  Phenomena  such  as  screening  of  the  carrier-phonon 
interaction,  carrier-carrier  scattering  are  thought  to  dominate  the  dissipation 
process  at  high  carrier  concentration.  Further  complications  might  arise  also 
as  a  result  of  plasmon-phonon  coupling.  It  is  clem  that  am  exact  theoretical 
treatment  of  such  a  system  is  today  out  of  range,  althought  a  big  effort  is  spent 
in  this  direction.  On  the  other  side,  experimental  measurements  usually  give  an 
integrated  view  that  might  not  allow  to  isolate  the  effects  of  specific  phenom¬ 
ena.  With  this  in  mind,  the  present  chapter  will  outline  in  Sec.  4.2  some  of  the 
theoretical  approaches  that  have  been  pursued,  together  with  the  details  of  the 
Monte  Carlo  algorithm.  The  application  of  the  Monte  Carlo  procedure  to  the 
case  of  laser  excitation  in  bulk  GaAs  and  GaAs/AlGaAs  quantum  wells  will  be 
discussed  in  Sec.  4.3.  Throughout  the  chapter,  the  Monte  Carlo  analysis  will 
only  be  limited  to  situations  in  which  the  adopted  physical  model  is  justified. 

4.2  The  Transport  Model  and  the  Monte  Carlo  algorithm 

The  dynamical  evolution  of  the  carrier  phonon  system  can  be  adeguately 
described  by  the  coupled  Boltzmann  equations: 

4- SI  +SI  +SI  +il  «■» 

dt  oilfield  at  Ic-ph  dt  lc— c  dt  \c— impurity 
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(4.2) 


dNq^dNq  |  dATq  I 

dt  dt  Ipk— e  dt  Ipk— ph 
where  /  and  JVq  are  respectively  the  carrier  and  the  phonon  distribution  func¬ 
tions.  In  the  following,  only  electrons  will  be  considered.  Such  an  assumption 
is  justified  as  long  as  the  carrier  densities  are  not  too  high,  so  that  no  effi¬ 
cient  coupling  exist  between  electrons  and  holes  (electron  concentrations  always 
lower  than  5  x  1017em~ 3  will  be  dealt  with  in  the  paper).  In  the  same  spirit, 
unscreened  electron-phonon  interaction  will  be  treated.  The  tim^-dependent 
transport  equations  for  carriers  and  phonons  are  coupled  through  the  occurence 
in  the  carrier-phonon  collision  integrals  of  both  carrier  and  phonon  distribution 
functions  [4.1).  A  decisive  simplification  of  the  phonon  equation  comes  from  the 
possibility  to  use  a  relaxation  time  for  the  phonon-phonon  interactions,  in  the 
form: 

dN{ q)  |  N (q)  —  ATj,  ^ 

dt  lph-ph  T0 p 

where  is  Ni  is  the  thermal  Plank  distribution 

NL  =  ('*"k'T‘-l)~1  (4.4) 

The  relaxation  time  approximation  is  justified  by  the  fact  that  the  phonon- 
phonon  interactions  are  dominated  by  the  decay  of  the  LO  phonons  into  pairs 
of  electronically  non  active  phonons  from  zone-boundary  modes,  values  of  the 
phonon  lifetime  rop  are  generally  of  the  order  of  10  ps,  with  a  weak  temperature 
dependence  [4.2). 

Time-resolved  phonon  spectroscopy  has  yield  a  rather  wide  range  of  values 
for  r,  between  7  ps  [4.3]  and  28  ps  [4.4].  The  reason  for  this  spread  of  experi¬ 
mentally  determined  LO  phonon  lifetimes  seems  to  have  two  sources.  Firstly, 
the  quality  of  the  sample  surface  can  strongly  influence  the  decay  dynamics 
within  the  thin  light-absorption  layer  [4.5].  Secondly,  the  decay  rate  of  a  non 
thermal  phonon  population  might  contain  strong  contributions  from  the  reab¬ 
sorption  by  the  photogenerated  carriers  of  the  Initially  excited  phonons.  This 


-68- 


point  will  be  discussed  in  detail  later.  Our  choice  of  r  equal  to  7  ps  at  77  K  and 
3.5  at  300  K  is  in  agreement  with  the  most  recent  experimental  results  [4.6,4.7]. 
Several  theoretical  approaches  have  been  presented  in  the  literature  for  the  so¬ 
lution  of  Eqs.  4.1  and  4.2.  Details  about  the  various  methods  can  be  found 
in  the  references.  Certainly  one  of  the  most  interesting  one  is  due  to  Collet 
and  Amand  [4.8]  who  directly  solved  the  coupled  transport  equations  through 
a  discretization  in  q-space  and  in  time  to  obtain  the  evolution  of  the  carrier 
and  phonon  distributions  and  of  the  mean  e-h  plasma  energy  during  and  after 
80  femtosecond  laser  excitation  pulses  of  varying  intensity.  Another  method, 
called  the  heated  and  drifted  Maxwellian  (HDM)  model,  is  of  particular  interest 
here  since  it  originated  the  Monte  Carlo  investigation  of  phonon  perturbation. 
In  the  HDM  approach,  the  carriers  are  assumed  to  be  characterized  by  a  heated 
and  drifted  Maxwellian  distribution.  In  the  first  attempt  to  use  a  Monte  Carlo 
technique  in  the  study  of  phonon  perturbations,  it  was  assumed  that  the  per¬ 
turbed  distribution  function  for  acoustic  phonons  of  wavevector  q  was  given  by 
the  asymptotic  steady-state  solution  N q  of  the  phonon  Boltzmann  equation  in 
the  presence  of  a  HDM  carrier  distribution 


N  =  ^c9c  + 

**  Se  +  Sph 


(4.5) 


where 

<«> 

a  Planck  distribution,  heated  to  a  carrier  temperature  Tc  and  shifted  about 
Vi.  In  Eq.  (4.5),  ge  and  yph  are  respectively  the  inverse  relaxation  time  for 
electronic  and  non  electronic  phonon  transitions.  At  a  given  temperature  N, q 
depends  basically  on  three  physical  parameters  which  describe  the  ensemble  of 
the  carriers,  i.e.:  the  concentration  ne,  the  mean  drift  velocity  vj  and  the  tem¬ 
perature  Tc  of  the  carriers.  These  parameters  contribute  in  different  ways  to  the 


determination  of  jVq.  An  increase  in  nc  will  lead  to  an  increase  of  the  phonon 
perturbation.  Moreover  an  increase  of  vs  will  increase  the  phonon  perturbation 
for  modes  with  wavectors  within  the  forward-cone  around  vj.  An  increase  of 
Tc  results  in  an  overall  increase  of  #q.  Finally,  lower  lattice  temperatures,  via 
Eq.(4.5),  also  favour  the  phonon  perturbation.  To  evaluate  the  effects  of  the 
phonon  perturbation  on  the  semiconductor  transport  properties,  an  iterative 
scheme  was  used.  Starting  a  MC  simulation  with  the  equilibrium  phonon  dis¬ 
tribution,  the  values  of  the  carrier  mean  energy  <  e  >  and  drift  velocity  Vd 
were  determined  and  substitute  in  Eq.  (4.1)  to  obtain  the  perturbed  phonon 
distribution.  The  new  scattering  rates,  accounting  for  the  phonon  perturbation, 
were  then  calculated  numerically  and  the  MC  algorithm  repeated  until  conver¬ 
gence  was  achieved.  A  semiconductor  model  parameterized  for  the  case  of  holes 
in  Ge  was  used.  The  effect  of  the  phonon  disturbance  leads  to  an  increase  of 
both  Vd  and  <  e  >  at  low  fields  and  to  a  decrease  of  these  quantities  at  higher 
fields.  The  reason  of  this  behavior  is  associated  with  the  competing  roles  of  the 
anisotropic  and  of  the  isotropic  contribution  to  the  phonon  amplification.  When 
the  anisotropic  aspect  prevails,  the  phonon-gas  drags  the  carriers  whose  drift 
velocity  and  mean  energy  increase  above  their  unperturbed  value  (directional 
effect  of  the  phonons).  When,  on  the  contrary,  the  net  isotropic  amplification 
dominates,  the  phonon  gas  most  effectively  randomizes  the  motion  of  the  car¬ 
riers;  therefore  Vd  and  <  e  >  are  reduced  with  respect  to  their  unperturbed 
values. 

Fig.  4.2  shows  the  drift  velocity  and  mean  energy  as  functions  of  the 
applied  electric  field  at  4.2  K  with  and  without  phonon  perturbation. 

The  method  just  described  was  still  lacking  consistency,  since  it  relied 
on  the  assumption  of  a  HDM  electron  distribution.  A  full  consistent  Monte 
Carlo  simulation  was  developed  for  the  case  of  photoexcitation  in  GaAs.  The 
novel  procedure  allows  to  follow  the  time  evolution  of  the  phonon  distribution. 
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Preliminary  results  have  been  presented  in  [4.9]. 

A  two-valley  (T  and  L  )  model  is  used  for  GaAs  (under  the  conditions  con¬ 
sidered  here  X  valleys  do  not  contribute  significantly).  The  following  scattering 
mechanisms  are  considered: 

-  acoustic  phonons  with  deformation  potential  coupling  ( Dac  —  7eV),  treated 
exactly  according  to  the  procedure  given  in  [4.10]; 

-  polar  optical  phonons,  without  screening; 

•  ionized  imputities,  treated  in  the  Conwell-Weisskopf  formalism; 

-  intervalley  T  —*  L  phonons  =  8x10 aeV/em);  Throughout  the  paper, 
nominally  undoped  materials  (with  a  residual  impurity  concentration  of  10M 
cm  3)  will  be  considered. 

-  electron-electron  interaction  between  T- valley  electrons  have  been  included 
using  the  algorithm  presented  in  Ref.  4.11.  There,  it  was  also  shown  that  at 
the  low  injection  densities  examined  in  the  present  work  (typically  5xl016  cm  3) 
,  the  effect  of  the  carrier-carrier  scattering  is  negligible. 

The  laser  excitation  is  reproduced  by  adding  particles  to  the  simulation, 
distributed  in  time  as  the  lineshape  of  the  laser  pulse,  as  shown  in  the  insert  of 
Fig.  4.3.  The  simulation  is  subdivided  in  time  intervals  At  (with  At  typically 
much  shorter  than  the  average  scattering  time  for  the  LO  phonon  scattering). 
At  time  T  =  /At,  the  number  of  Monte  Carlo  electrons  is  updated  from  the 
previous  step  according  to  the  espression  : 

N(T)  =  N(T  -  At)  +  C  At  cosh~l(uT),  (4.7) 

where  w  and  C  are  parameters  related  to  the  width  and  power  of  the  laser  pulse. 

Electrons  are  excited  in  the  conduction  band  centered  around  a  given  en¬ 
ergy  Einj  ,  with  a  small  broadening  depending  on  the  width  of  the  laser  pulse. 
Since  the  excitation  energies  considered  here  are  below  the  threshold  for  in¬ 
tervalley  scattering  (0.3  eV  for  T  to  L  transitions),  there  is  here  no  significant 
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transfer  to  the  satellite  valleys. 

The  LO-phonon  distribution  function  is  followed  in  its  time  evolution, 
and  phonon-induced  modifications  to  the  relaxation  rates  of  the  electrons  are 
considered.  The  disturbances  of  other  types  of  phonons  are  negligible  in  the 
situation  examined  here.  In  a  finite  difference  scheme,  Eq.  (4.2)  for  the  phonon 
evolution  can  be  written  in  the  form: 

jV,(nAi)  =  AT,((n  -  1)  At)  +  dtf,(nAt)|  _ 

,  x  At  (4-8) 

-  (Nq(nAt)  -  N0) —  ;  n  =  1,2,3... 

The  procedure  set  up  to  account  for  the  LO-phonon  disturbances  has  the  fol¬ 
lowing  features: 

i)  the  time  evolution  of  the  LO-phonon  distribution  JVq  is  calculated  as  a  func¬ 
tion  of  wavevector  q  from  the  Monte  Carlo  simulation  ,  by  setting  up  a  his¬ 
togram  hq  defined  over  a  grid  in  q-space  of  mesh  size  A q.  After  each  scattering 
event  involving  an  LO  phonon,  the  histogram  is  updated.  In  the  absence  of 
external  d.c.  fields,  because  of  the  full  spherical  symmetry  only  the  amplitude 
of  q  is  relevant,  thus  reducing  the  complexity  and  the  storage  requirements 
of  the  simulation.  Preliminary  results  for  the  field  dependent  case  have  been 
presented  recently  [4.12]. 

ii)  At  fixed  times  T  =  j At  during  the  simulation  ,  Nq  is  calculated  as 

Nq{jAt)  =  Nq(jAt)  +  [Nq(jAt)  -  JVbl^.  (4.9) 

Here,  No  is  the  equilibrium  Bose  distribution  and 

ffq{}  At)  =  Nq({]  -  1)  At)  +  AAhq .  (4.10) 

The  term  A  A hq  gives  the  dynamical  contribution  of  the  electronic  processes  to 
the  phonon  distribution  during  the  time  step  At.  There,  A  is  a  normalization 


factor  accounting  for  the  density  of  states  in  q  -space  and  for  the  concentration 
of  excited  electrons,  given  by 


2r  n 


(4.11) 


where  n  is  the  electron  concentration  and  N  the  number  of  simulated  particles. 
The  second  term  on  the  right  hand  side  of  Eq.  (4.9)  accounts  for  the  phonon- 
phonon  processes.  The  algorithm  for  the  phonon  counting  can  be  viewed  as  a 
hybrid  Monte  Carlo  solution  of  the  phonon  Boltzmann  equation  within  a  finite 
difference  scheme. 

iii)  to  account  for  the  modifications  induced  by  the  phonon  disturbance  on  the 
rate  of  electron-phonon  scatterings,  the  integrated  scattering  probabilities  for 
LO-phonons  are  calculated  and  tabulated  at  the  beginning  of  the  simulation 
using  an  artificially  high  value  Nmax  for  the  phonon  distribution.  The  choice  of 
the  final  state  of  each  scattering  process  involving  an  LO  phonon  is  made  using  a 
rejection  technique  which  compares  the  actual  value  of  the  differential  scattering 
rate  with  the  maximized  one.  In  this  way,  we  are  able  to  discriminate  between 
the  scatterings  that  can  be  attributed  to  the  enhanced  phonon  distribution 
versus  those  induced  by  the  initial  maximization  (  which  axe  treated  as  self- 
scatterings  in  the  simulation) .  In  order  to  reduce  the  number  of  self  scattering 
events  ,  it  is  possible  to  recalculate  the  scattering  rates  at  fixed  times  during 
the  simulation.  A  numerical  integration  over  the  perturbed  phonon  distribution 
function  gives  the  exact  scattering  rates  at  a  given  time,  which  can  be  used 
directly  on  the  simulation. 

The  suggested  procedure  is  a  full  Monte  Carlo  simulation  of  the  dynamics 
of  an  interacting  electron-phonon  gas  within  the  finite  difference  scheme  indi¬ 
cated  above,  free  of  adjustable  parameters.  In  the  next  section,  the  results  of 
the  Monte  Carlo  simulation  in  the  presence  of  laser  excitation  will  be  presented. 
The  application  of  the  algorithm  in  the  precence  of  an  applied  electric  field  will 
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be  presented  in  a  separate  contribution.  In  that  case,  a  two  dimensional  grid  is 
used  to  describe  the  perturbed  diribution. 

4.3  Applications 

The  question  of  possible  effects  of  nonequilibrium  optical  phonon  distri¬ 
butions  on  the  dynamics  of  optically  excited  charge  carriers  in  semiconduc¬ 
tors  is  becoming  a  widely  investigated  and  debated  topic.  The  main  scientific 
motivation  came  from  the  rapid  development  of  picosecond  and  subpicosecond 
laser  spectroscopy,  which  allows  to  study  even  the  fastest  relaxation  phenomena 
in  solids  and  thereby  also  some  fundamental  hot  carrier-hot  phonon  processes 
which  might  ultimately  limit  the  switching  efficiencies  of  uitrafast  electronic 
devices.  The  algorithm  described  in  the  previous  section  has  been  applied  to 
various  situations  to  study  the  dynamics  of  the  LO  phonon,  of  the  electron 
distributions,  and  their  mutual  effects.  The  time  evolution  of  the  perturbed 
phonon  distribution  is  shown  in  Fig.  4.3  for  an  excited  carrier  density  of  5xl016 
cm-*.  Electrons  are  excited  at  0.25  eV  above  the  bottom  of  the  conduction 
band,  corresponding  to  a  photon  energy  of  1.8  eV.  The  lattice  temperature  is 
77  K.  The  lineshape  of  the  laser  pulse  is  shown  in  the  insert  (halfwidth  =  0.8 
ps).  The  LO  distribution  is  driven  out  of  equilibrium  even  during  the  excitation, 
due  to  the  fast  power  dissipation  of  the  high-energy  photoexcited  electrons.  The 
maximum  is  reached  at  a  delay  time  of  1  picosecond  for  wavevectors  of  about 
6x10*  cm  '.  The  small  q  values  that  are  amplified  during  and  immediately 
after  the  excitation  are  due  to  the  polar  nature  of  the  e-phonon  coupling.  At 
longer  times,  the  phonon  distribution  relaxes  towards  its  equilibrium  value  as  a 
result  of  two  distinct  processes,  phonon  reabsorption  and  phonon-phonon  inter¬ 
action.  The  first  one  is  due  to  the  fact  that  the  group  velocity  of  optical  phonons 
is  very  small  (less  than  103cm/a),  implying  that  the  phonons  cannot  drift  away 
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of  the  excitation  volume  during  their  lifetime.  Therefore,  if  the  phonon  lifetime 
is  long  enough  and  the  coupling  with  the  carriers  sufficiently  strong,  emitted 
phonons  can  be  reabsorbed. 

It  is  important  to  notice  that  modes  of  different  wavevector  evolve  in 
time  in  different  ways,  as  indicated  in  Fig.  4.4  .  Those  with  the  smaller  q 
(6x10s  cm-1)  exhibit  an  exponential  decay,  immediately  after  the  end  of  the 
excitation,  with  a  characteristic  decay  time  of  7  ps.  At  intermediate  q’s  (8  and 
10x10s  em~  *)  the  phonon  distribution  decay  much  faster  at  short  times  (up  to  5 
and  8  ps  delay),  approaching  then  the  exponential  behavior.  The  amplification 
of  these  large-q  phonons  is  not  as  pronounced  as  that  of  the  small-q  ones. 

The  time  evolution  of  the  phonon  distribution  reflects  the  microscopic  de¬ 
tails  of  the  cooling  processes  in  the  coupled  electron-phonon  system.  While 
phonon-phonon  processes  are  always  active,  and  their  effect  is  independent  of 
wavevector,  phonon  reabsorption  varies  drastically  as  a  function  of  time  and 
wavevector.  In  fact,  the  very  rapid  changes  in  the  electron  distribution  func¬ 
tion  (that  will  be  examine  below)  modify  the  range  of  phonon  transitions  that 
are  allowed  by  energy  and  momentum  conservation.  Fig.  4.5  shows  the  mini¬ 
mum  q  for  LO-phonon  absorption  and  emission  as  a  function  of  electron  energy 
in  a  parabolic  band.  At  high  energy,  electrons  can  emit  phonons  with  very  small 
q,  but  as  they  cool  the  minimum  allowed  q  shifts  at  higher  values.  Such  a  shift 
appears  in  Fig.  4.3,  hidden  though  by  the  strong  initial  amplification.  Further¬ 
more,  an  electron  will  not  be  able  to  reabsorb  the  earlier  emitted  phonons  once 
it  goes  below  a  certain  energy. 

This  simple  analysis  explains  why  the  phonons  with  small  q- vector  excited 
during  the  first  stages  of  the  electron  relaxation  (up  to  2  ps  delay  time)  cannot 
be  reabsorbed,  and  decay  exponentially  via  non-electronic  processes. 

On  the  other  side,  both  the  reabsorption  and  the  phonon- phonon  terms 
will  contribute  to  the  damping  of  phonons  of  larger  wavevector  in  the  first  few 
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picoseconds,  leading  to  their  faster  decay  over  this  time  interval. 

The  modification  of  the  scattering  rates  for  the  electron-LO  phonon  in¬ 
teraction  due  to  the  phonon  perturbation  is  presented  in  Fig.  4.6.  There  the 
total  scattering  rates  for  absorption  and  emission,  obtained  from  a  numerical 
integration  over  the  perturbed  distribution  function  Nq,  are  plotted  at  different 
time  delays  after  the  end  of  the  laser  pulse.  Since  at  low  temperature  the  value 
of  the  equilibrium  phonon  distribution  is  much  smaller  than  unity,  the  emission 
probability  is  a  few  orders  of  magnitude  higher  than  the  absorption  one.  As 
the  phonon  population  grows  out  of  equilibrium,  the  absorption  rate  increases 
dramatically,  relatively  much  faster  than  the  emission  one.  The  changes  of  the 
scattering  rate  with  time  reflect  the  temporal  evolution  of  the  phonon  popu¬ 
lation.  It  is  important  to  notice  that  even  a  few  picoseconds  after  the  pulse, 
a  significant  amount  of  phonons  is  still  present  and  a  considerable  number  of 
phonon  reabsorptions  are  detected. 

The  time  evolution  of  the  electron  distribution  function,  shown  in  Fig. 
4.7,  completes  the  previous  analysis  of  the  phonon  amplification.  The  distinct 
peaks  in  the  distributions  at  short  time  delays  (0  and  1  ps)  are  due  mainly  to 
LO  phonon  emission  which  sets  up  already  during  the  laser  pulse  (an  average 
time  of  160  fs  for  the  emission  of  an  LO  phonon  by  electrons  at  the  excitation 
energy  is  calculated  from  the  simulation).  At  a  time  delay  of  4  ps,  the  electrons 
mainly  populate  the  low  energy  region  below  100  meV.  Many  of  them  have  an 
energy  below  the  threshold  for  LO  phonon  emission.  It  will  be  seen  later  that 
in  this  case  reabsorption  can  become  very  important. 

As  a  last  remark  on  the  phonon  dynamics,  it  is  important  to  compare  the 
previous  considerations  with  the  experimental  results  of  Raman  spectroscopy. 
The  shaded  area  in  Fig.  4.5  indicates  the  range  of  raman-active  wavevectors  for 
the  data  given  in  Refs.  (4.6,4.13].The  Monte  Carlo  result  for  those  modes  (curve 
x  in  Fig.  4.4)  are  in  good  qualitative  agreement  with  the  findings  of  Hash  et 
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al.  [4.6],  obtained  for  the  same  electron  concentration  but  a  higher  excitation 
energy. 

The  effect  of  the  phonon  perturbation  on  the  cooling  of  the  photoexcited 
electrons  is  shown  in  Fig.  4.8.  The  electron  relaxation  rate  is  drastically  reduced 
because  of  the  presence  of  non-equilibrium  phonons.  The  phenomenon  is  mainly 
due  to  the  reabsorption  of  the  LO  phonons  that  have  been  emitted  in  the  first 
stage  of  the  relaxation  without  having  had  enough  time  to  decay.  The  effect  of 
phonon  reabsorption  grows  with  time  as  the  electrons  populate  the  low  energy 
regions  below  the  threshold  for  optical  emission.  It  has  also  been  found  that  the 
reduction  in  the  cooling  rate  of  the  electrons  is  even  larger  at  higher  electron 
densities  or  higher  injection  levels. 

The  Monte  Carlo  algorithm  presented  here  has  been  compared  with  the  models 
of  Ref.  4.14  and  4.5.  At  the  low  excitation  energies  used  in  Ref.  4.15,  the 
phonon  disturbance  is  reduced  with  respect  to  the  case  shown  if  Fig.  4.3.,  and 
reaches  its  maximum  at  higher  q’s.  The  Monte  Carlo  result  agrees  quite  well 
with  ones  of  the  more  sophisticated  model  of  Collet  et  al. 

The  temperature  model  of  Ref.  4.14  assumes  that  the  carriers  (electrons 
and  holes)  are  characterized  by  a  Fermi-Dirac  distribution  at  any  time  during 
and  after  the  laser  pulse,  corresponding  to  a  very  fast  thermalization  within 
the  photogenerated  plasma.  In  order  to  verify  the  consistency  of  our  results, 
we  have  performed  a  simulation  by  assuming  that  the  carriers  are  initially  dis¬ 
tributed  according  to  a  heated  maxwellian  distribution.  The  Monte  Carlo  re¬ 
sults  obtained  using  the  same  parameters  of  Fig.  4.3  indicate  in  this  case  a 
much  smaller  phonon  perturbation,  with  the  maximum  of  the  phonon  distri¬ 
bution  still  reached  at  1  picosecond  delay  time  as  in  Fig.  4.3,  but  its  value 
is  reduced  by  a  factor  two.  The  reduction  of  phonon  heating  is  related  to  the 
increased  population  of  the  low  energy  tails  of  the  maxwellian  distribution  com¬ 
pared  to  the  quasi  monoenergetic  distribution  used  before.  The  Monte  Carlo 
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results  agree  very  well  with  those  of  the  temperature  model,  for  a  one  compo¬ 
nent  system.  When  both  electrons  and  holes  are  considered,  the  temperature 
model  (which  assumes  the  same  temperature  for  both  carriers)  shows  a  further 
reduction  of  phonon  heating,  which  indicates  a  strong  energy  transfer  from  the 
electron  to  the  hole  system,  which  then  dissipate  mainly  via  TO  emission.  Al¬ 
though  the  latter  result  depends  heavily  on  the  assumptions  of  the  model,  it 
nevertheless  shows  that  the  e-h  interaction  can  be  very  important.  A  prelimi¬ 
nary  step  to  combine  the  effect  of  e-h  scattering  and  non  equilibrium  phonons 
has  been  recently  presented  [4.16].  The  results  just  shown  refers  to  the  simple 
case  of  GaAs,  where  only  the  central  valley  is  important.  In  general,  especially 
if  the  excitation  energy  is  sufficiently  high,  the  population  of  the  higher  val¬ 
leys  (L  and  X)  is  not  negligible.  The  influence  of  such  effect  on  the  phonon 
population  can  be  very  strong.  Fig.  4.9  shows  the  minimum  q  for  LO  phonon 
emission  as  a  function  of  electron  energy  for  the  T  (same  curve  as  Fig.  4.4) 
and  the  L  valleys.  Due  to  the  higher  effective  mass  of  the  satelite  valley,  the 
emitted  LO  phonons  have  a  larger  wavevector.  Since  the  area  of  phase-space  is 
increased,  their  contribution  to  the  phonon  population  will  be  redufed.  This  is 
illustrated  in  Fig.  4.10,  where  an  initial  electron  energy  of  0.5  eV  has  been  used 
(with  a  T  —  L  separation  of  0.3  eV).  The  parameters  of  the  simulation  are  the 
same  as  before.  After  l  picosecond  from  the  excitation,  about  60  percent  of  the 
electrons  are  found  in  the  L  valley.  The  Monte  Carlo  hystogram  (Fig.  4.10a) 
confirms  that  the  emission  of  LO  phonons  by  L*vaIIey  electrons  is  concentrated 
in  the  large  q  region.  The  actual  number  of  phonons  reflects  indeed  the  rel¬ 
ative  population  of  the  two  valleys.  Nevertheless,  the  effect  of  those  phonons 
on  the  perturbed  distibution  (Fig.  4.10b)  is  negligible.  Furthermore,  all  of  the 
L-valley  phonon  have  q  values  too  large  to  be  detected  spectroscopically.  We 
can  therefore  expect  that  phonon  amplifications  experimentally  detected  would 
decrease  when  a  relevant  number  of  intervalley  transfers  is  present. 
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Recent  results  obtained  with  time  resolved  photoluminescence  have  shown 
dramatically  reduced  energy  relaxation  rates  for  photo-excited  electrons  in 
GaAs-AlGaAs  quantum  wells  [4.17 , 4.18j.  It  is  still  an  open  question  whether 
the  slow  cooling  rates  are  due  to  the  presence  of  non-equilibrium  phonons  or  to 
the  effect  of  reduced  dimensionality  and  screening  [4.19,4-20].  The  algorithm 
just  presented  has  been  applied  to  a  single  quantum  well  of  GaAs-AlGaAs  (tw 
=  150  A,  4>b  =  0.28  eV),  with  subband  energies  given  by  the  solution  of  the 
one  dimensional  wave  equation  for  a  square  well  potential.  The  bands  are  as¬ 
sumed  parabolic.  The  scattering  rates  of  the  quantized  2D  electrons  with  bulk 
unscreened  LO-phonon  (both  intra-  and  intersubband)  are  calculated  numeri¬ 
cally  without  the  use  of  momentum  conserving  approximations  [4.21].  It  has 
been  shown  [4.21]  that,  for  wells  larger  than  100  A,  very  little  difference  exist 
from  the  scattering  rates  calculated  accounting  for  phonon  confinement  (slab 
modes)  and  the  one  obtained  using  bulk  modes.  Intervalley  transfer  to  the 
L-valleys  (also  quantized)  is  included  as  well.  2D  electron-electron  scattering 
is  introduced  in  the  Monte  Carlo  simulation  through  a  generalization  of  the 
self-scattering  technique  given  in  Ref.  4.22  to  the  multi-subband  quantized 
system  [4.23].  The  various  electrons  are  allowed  to  interact  via  a  statically 
screened  Coulomb  interaction  determined  by  the  long  wavelength  limit  of  the 
twodimensional  Lindhard  dielectric  function.  Degeneracy  effects  due  to  the 
Pauli  exclusion  principle  are  also  considered  [4.24].  In  the  present  simulation 
we  have  neglected  electron-hole  scattering  and  recombination,  which  might  be 
of  importance  in  some  of  the  reported  experiments. 

In  order  to  include  non-equilibrium  LO  phonons  in  the  Monte  Carlo  sim¬ 
ulation,  we  generalize  to  twodimensional  systems  the  procedure  described  ear¬ 
lier  for  bulk  GaAs.  The  phonon  distribution  is  given  directly  by  a  detailed 
balance  of  the  emission  and  absorption  events  during  the  simulation  on  a  grid 
in  q-space  with  the  excess  LO  phonon  population  in  each  mode  decaying  via  a 
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phenomenological  phonon  lifetime  of  7  ps.  Since  no  direct  measurement  of  LO 
phonon  lifetime  in  quantum  well  systems  is  available  at  this  time,  the  measured 
bulk  value  is  used.  In  2D,  the  component  of  the  phonon  wavevector  in  the  di¬ 
rection  of  the  well,  qt,  is  not  conserved.  For  a  square  well,  however,  it  has  been 
shown  that  the  qz  component  is  very  peaked  for  wavevectors  corresponding  to 
the  change  in  intersubband  energy  during  the  transition  [4.14].  Therefore,  the 
phonon  distribution  is  tabulated  for  discrete  qz  corresponding  to  the  various 
inter-  and  intrasubband  events,  with  the  phonon  wavevectors  otherwise  treated 
as  two-dimensional. 

'  The  cooling  of  photoexcited  electrons  in  an  n-type  GaAs-AlGaAs  quantum 
well  at  low  temperature  (5  K)  has  been  considered  [4.25,  4.26].  A  background 
density  of  2.5xlOu  cm-2  is  used.  The  injected  density  is  5xlOn  cm-2  The 
GaAs  parameters  are  the  same  as  for  the  bulk  case.  The  width  of  the  simulated 
laser  pulse  is  about  1  ps,  during  which  time  carriers  are  added  to  the  simulation 
with  an  initial  energy  of  0.25  eV  above  the  bottom  of  the  lowest  subband.  Fig. 
4.11  shows  the  evolution  of  the  electron  total  energy  (kinetic  plus  potential)  as 
a  function  of  time  during  and  after  the  pulse.  The  excited  electrons  lose  energy 
mainly  through  the  interaction  with  the  background  electrons  and  through  the 
emission  of  LO  phonons.  Without  hot  phonons,  the  hot  electrons  are  found  to 
reach  equilibrium  in  about  3  ps.  In  contrast  with  the  case  of  an  unperturbed 
phonon  distribution,  a  much  slower  relaxation  is  found  when  non-equilibrium 
phonons  are  accounted  for.  The  two  cases  in  Fig.  4.11  correspond  to  different 
quantum  wells.  It  is  clear  that  no  real  dependence  on  well  width  has  been  found 
in  the  Monte  Carlo  result.  The  total  electron  energy  plotted  here  is  calculated  as 
an  ensemble  average  during  the  simulation.  For  degenerate  systems,  this  quntity 
can  vary  considerably  from  the  electron  temperature,  which  can  be  rigorously 
defined  only  in  the  presence  of  a  fermi  distribution,  the  elctron  temperature 
calculated  from  the  average  energy  ,  assuming  the  distribution  function  is  fermi- 
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like,  (dashed  curve)  correspond  to  the  slope  of  the  tail  of  the  distribution. 
The  reduction  of  the  electron  cooling  rate  is  due  to  the  reabsorption  of  non¬ 
equilibrium  phonons  which  build  up  during  the  initial  pumping  and  the  first 
stage  of  the  electron  relaxation.  This  is  evidenced  by  the  time  evolution  of  the 
phonon  distribution  at  qz= 0  (intra  subband  scattering)  shown  in  Fig.  4.12.  The 
same,  although  reduced,  features  are  found  also  at  qz= 0.  LO-phonon  emission 
during  the  pulse  and  immediately  after  creates  a  large  population  of  phonons 
at  small  q’s.  At  longer  times,  phonon  reabsorptions  and  phonon-phonon  losses 
drive  the  distribution  back  to  equilibrium.  The  secondary  peak  that  develops 
at  later  times  in  the  phonon  distribution  is  due  to  phonon  emission  by  electrons 
that  have  already  relaxed  to  lower  energy. 

As  pointed  out  before,  the  reduction  in  the  electron  relaxation  rate  is 
mainly  due  to  reabsorption  of  the  emitted  LO-phonons.  The  effect  is  stronger 
when  a  considerable  number  of  electrons  have  relaxed  to  the  low  energy  region 
below  the  emission  threshold. 

Shortly  after  the  end  of  the  laser  pulse  (t  =  1.6  ps),  the  strong  intercarrier 
scattering  creates  a  broad  distribution  where  the  subband  minima  (indicated  by 
arrows)  clearly  appear.  Within  each  subband,  the  distribution  function  starts 
to  exhibit  a  Fermi-like  appearance  wich  is  fully  established  at  longer  times  as 
shown  in  the  bottom  curve. 

4.4  Conclusions 

In  summary,  a  full  Monte  Carlo  technique  for  the  study  of  electron  and 
phonon  dynamics  in  bulk  GaAs  and  GaAs-AlGaAs  quantum  wells  has  been 
presented.  A  strongly  perturbed  phonon  distribution  is  found  in  the  first  pi¬ 
coseconds  after  the  laser  pulse,  which  is  responsible  for  a  reduction  of  the  cooling 
rate  of  the  photoexited  carriers.  Good  qualitative  agreement  with  Raman  data 
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Figure  Captions 


Fig.  4.1  Schematic  representation  of  the  energy  flux  inside  and  outside  a  semi¬ 
conductor. 

Fig.  4.2  Drift  velocity  and  mean  energy  as  functions  of  the  electric  field  strength 
for  holes  in  Gc  at  4.2  K  for  the  indicated  impurity  concentration  Ni.  Dashed 
(continuous)  carves  refer  to  an  unperturbed  (perturbed)  phonon  distribution 
function. 

Fig.  4.3  Non-equilibrium  LO  phonon  distribution  functions  at  three  different 
delay  times  as  afunction  of  the  phonon  wavevector.  The  insert  shows  the  shape 
of  the  laser  pulse. 

Fig.  4.4  Time  evolution  of  four  different  modes  as  a  function  of  delay  time. 

Fig.  4.5  Energy  dependence  of  the  minimum  wavevector  for  absorption  and 
emission  of  LO  phonons,  as  a  function  of  energy  in  a  parabolic  band.  The 
shaded  are  indicates  typical  raman  active  phonon  wavevectors. 

Fig.  4.6  Absorption  and  emission  scattering  rates  for  LO  phonons  as  a  function 
of  energy,  at  different  times  after  the  laser  pulse. 

Fig.  4.7  Time  evolution  of  the  electron  distribution  function,  at  different  delay 
times. 

Fig.  4.8  Average  electron  energy  as  a  function  of  time  with  (continuous  curve) 
and  without  (dashed  curve)  hot-phonons. 

Fig.  4.9  Energy  dependence  of  the  minimum  wavevector  for  absorption  and 
emission  of  LO  phonons,  as  a  function  of  energy  in  a  parabolic  hand  for  T  and 
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